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Recent results on theoretical studies of heat conduction in low-dimensional systems are pre- 
sented. These studies are on simple, yet nontrivial, models. Most of these are classical systems, 
but some quantum-mechanical work is also reported. Much of the work has been on lattice 
models corresponding to phononic systems, and some on hard particle and hard disc systems. 
A recently developed approach, using generalized Langevin equations and phonon Green's 
functions, is explained and several applications to harmonic systems are given. For interact- 
ing systems, various analytic approaches based on the Green-Kubo formula are described, and 
their predictions are compared with the latest results from simulation. These results indicate 
that for momentum-conserving systems, transport is anomalous in one and two dimensions, 
and the thermal conductivity k, diverges with system size L, as k ~ L a . For one dimensional 
interacting systems there is strong numerical evidence for a universal exponent a = 1/3, but 
there is no exact proof for this so far. A brief discussion of some of the experiments on heat 
conduction in nanowires and nanotubes is also given. 
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1. Introduction 

It is now about two hundred years since Fourier first proposed the law of heat 
conduction that goes by his name. Consider a macroscopic system subjected to 
different temperatures at its boundaries. One assumes that it is possible to have a 



* Corresponding author. Email: dabhi@rri.res.in 



November 19, 2008 19:56 Advances in Physics reva 
2 

coarse-grained description with a clear separation between microscopic and macro- 
scopic scales. If this is achieved, it is then possible to define, at any spatial point x 
in the system and at time t, a local temperature field T(x, t) which varies slowly 
both in space and time (compared to microscopic scales). One then expects heat 
currents to flow inside the system and Fourier argued that the local heat current 
density J(x, t) is given by 

J(x,t) = -reVT(x,t) , (1) 

where k is the thermal conductivity of the system. If u(x, t) represents the local 
energy density then this satisfies the continuity equation du/dt + V.J = 0. Using 
the relation du/dT = c, where c is the specific heat per unit volume, leads to the 
heat diffusion equation: 

— L-l^ = ~V. «VT(x,t) 
at c 

Thus, Fourier's law implies diffusive transfer of energy. In terms of a microscopic 
picture, this can be understood in terms of the motion of the heat carriers, i.e. 
, molecules, electrons, lattice vibrations(phonons), etc., which suffer random col- 
lisions and hence move diffusively. Fourier's law is a phenomological law and has 
been enormously succesful in providing an accurate description of heat transport 
phenomena as observed in experimental systems. However there is no rigorous 
derivation of this law starting from a microscopic Hamiltonian description and 
this basic question has motivated a large number of studies on heat conduction in 
model systems. One important and somewhat surprising conclusion that emerges 
from these studies is that Fourier's law is probably not valid in one and two di- 
mensional systems, except when the system is attached to an external substrate 
potential. For three dimensional systems, one expects that Fourier's law is true in 
generic models, but it is not yet known as to what are the neccessary conditions. 

Since one is addressing a conceptual issue it makes sense to start by looking at 
the simplest models which incorporate the important features that one believes 
are necessary to see normal transport. For example, one expects that for a solid, 
anharmonicity and disorder play important roles in determining heat transport 
properties. Thus most of the theoretical studies have been on these simple models, 
rather than on detailed models including realistic interparticle potentials, etc. The 
hope is that the simple models capture the important physics, and understanding 
them in detail is the first step towards understanding more realistic models. This 
review almost exclusively will talk about simple models of heat conduction in low 
dimensional systems, mostly one dimensional (ID) and some two dimensional (2D). 
Also a lot of the models that have been studied are lattice models, where heat is 
transported by phonons, and are relevant for understanding heat conduction in 
electrically insulating materials. Some work on hard particle and hard disc systems 
will also be reviewed. 

There are two very good earlier review articles on this topic, including those by 
Bonetto et al. [l| and Lepri et al. Some areas that have not been covered 
in much detail here can be found in those reviews. Another good review, which 
also gives some historic perspective, is that by Jackson [3]. Apart from being an 
update on the older reviews, one area which has been covered extensively in this 
review is the use of the nonequilibrium Green's function approach for harmonic 
systems. This approach nicely shows the connection between results from various 
studies on heat transport in classical harmonic chain models, and results obtained 
from methods such as the Landauer formalism, which is widely used in mesoscopic 
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physics. As we will see, this is one of the few methods where explicit results can 
be obtained for the quantum case also. 

The article is organized as follows. In Sec. ((21), some basic definitions and a de- 
scription of some of the methods used in transport studies is given. In Sec. ([3]), 
results for the harmonic lattice are given. The nonequilibrium Green's function 
theory will be developed using the Langevin equation approach and various ap- 
plications of this method are described. The case of interacting particles ((non- 
harmonic inter-particle interactions) in one dimension is treated in Sec. ([4]). This 
section briefly summarizes the analytic approaches, and then gives results of the 
latest simulations in momentum-conserving and momentum non-conserving one 
dimensional systems. The next section [Sec. ©] looks at the joint effect of disorder 
and interactions in one dimensional systems. In Sec. © results on two dimensional 
interacting systems are presented while Sec. ([7]) gives results for billiard like systems 
of noninteracting particles. Some of the recent experimental results on nanowires 
and nanotubes are discussed in Sec. ©. Finally the conclusions of the review are 
summarized in Sec. ([9]) and a list of some interesting open problems is provided. 

2. Methods 

The most commonly used approaches in heat transport studies have been: (i) those 
which look at the nonequilibrium steady state obtained by connecting a system 
to reservoirs at different temperatures, and (ii) those based on the Green-Kubo 
relation between conductivity and equilibrium correlation functions. In this section 
we will introduce some of the definitions and concepts necessary in using these 
methods [sees. (|2.1)2.2p ]. Apart from these two methods, an approach that has 
been especially useful in understanding ballistic transport in mesoscopic systems, 
is the nonequilibrium Green's function method and we will describe this method in 
sec. (I2.3p . Ballistic transport of electrons refers to the case where electron-electron 
interactions are negligible. In the present context ballistic transport means that 
phonon-phonon interactions can be neglected. 

2.1. Heat bath models, definitions of current, temperature and conductivity 

To study steady state heat transport in a Hamiltonian system, one has to connect it 
to heat reservoirs. In this section we will first discuss some commonly used models 
of reservoirs, and give the definitions of heat current, temperature and thermal 
conductivity. It turns out that there are some subtle points involved here and we 
will try to explain these. 

First let us discuss a few models of heat baths that have been used in the litera- 
ture. For simplicity we here discuss the ID case since the generalization to higher 
dimension is straight-forward. We consider a classical ID system of particles in- 
teracting through a nearest neighbour interaction potential U and which are in an 
external potential V. The Hamiltonian is thus: 

N 2 N - 1 

H = E[5r + 1 + E u te ~ ( 3 ) 
i=i 1 i=i 

where {mi,xi,pi = mixi} for I = 1,2, ...N denotes the masses, positions and mo- 
menta of the N particles. For the moment we will assume that the interparticle 
potential is such that the particles do not cross each other and so their ordering 
on the line is maintained. 
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To drive a heat current in the above Hamiltonian system, one needs to connect 
it to heat reservoirs. Various models of baths have been used in the literature and 
here we discuss three popular ones. 

(i) Langevin baths: These are defined by adding additional force terms in the 
equation of motion of the particles in contact with baths. In the simplest form, 
the additional forces consist of a dissipative term, and a stochastic term, which 
is taken to be Gaussian white noise. Thus with Langevin reservoirs connected to 
particles I = 1 and I = N, the equations of motion are given by: 

Pi = h - —Pi + VL(t) 

mi 

pi = fl for Z = 2,3..JV-1 

Pn = In Pn + VR(t) ( 4 ) 

m N 

where fi = 

dxi 

is the usual Newtonian force on the I th particle. The noise terms given by r?L,_R are 
Gaussian, with zero mean, and related to the dissipation coefficients ^l,r by the 
usual fluctuation dissipation relations 

(VL(t) VL (t')) = 2k B T LlL 5(t - t') 
(VR(t)m(t')) = 2k B T RlR 5(t - t') 
(VL(t) m (t')) = , 

where T^,Tr are the temperatures of the left and right reservoirs respectively. 

More general Langevin baths where the noise is correlated will be described in 
sec. f)3.2f) . Here we briefly discuss one particular example of a correlated bath, 
namely the Rubin model. This model is obtained by connecting our system of 
interest to two reservoirs which are each described by semi-infinite harmonic oscil- 
lator chains with Hamiltonian of the form H b = XXi P i I 2 + YaLq( x i ~ x i+i) V 2 , 
where {Xi, P{\ denote reservoir degrees of freedom and Xq = 0. One assumes that 
the reservoirs are initially in thermal equilibrium at different temperatures and 
are then linearly coupled, at time t = — oo, to the two ends of the system. Let us 
assume the coupling of system with left reservoir to be of the form — x±Xi. Then, 
following the methods to be discussed in sec. (|3,2p . one finds that the effective 
equation of motion of the left-most particle is a generalized Langevin equation of 
the form: 

Pi = fi + f dt'E L (t - 0*1 (*') + VL(t) , (5) 

<J — CO 

where the fourier transform of the kernel is given by: 

t L {u) = / dt^ L {t)e iujt = e iq for \u\ < 2 
Jo 

= —e~ u for \lu\ > 2 , 



and q,v are defined through cos(g) = 1 — lo 2 /2, cosh(^) = uj 2 /2 — 1 respectively. 
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The noise correlations are now given by: 

(fj L (u)v L (u')) = — JmpiM] 8{u + J), (6) 

TTUJ 

where = (l/2ir) dtr/i(t)e w *. Similar equations of motion are obtained 

for the particle coupled to the right reservoir. 

(ii) Nose-Hoover baths: These are deterministic baths with time-reversible 
dynamics which however, surprisingly, have the ability to give rise to irreversible 
dissipative behaviour. In its simplest form, Nose-Hoover baths attached to the end 
particles of the system described by the Hamiltonian Eq. ([2D, are defined through 
the following equations of motion for the set of particles: 

Pi = h - Clpi 

pi = ft for Z = 2,3...iV-l 

Pn = In - CrPn , (7) 

where £l and £i? are also dynamical variables which satisfy the following equations 
of motion: 

£ R = J_ ( Pn _ ! 
Or \m N k B T R 

with Ol and Or as parameters which control the strength of coupling to reservoirs. 

Note that in both models (i) and (ii) of baths, we have described situations where 
baths are connected to particular particles and not located at fixed positions in 
space. These are particularly suited for simulations of lattice models, where par- 
ticles make small displacements about equilibrium positions. Of course one could 
modify the dynamics by saying that particles experience the bath forces (Langevin 
or Nose-Hoover type) whenever they are in a given region of space, and then these 
baths can be applied to fluids too. Another dynamics where the heat bath is located 
at a fixed position, and is particularly suitable for simulation of fluid systems, is 
the following: 

(iii) Maxwell baths: Here we take particles described by the Hamiltonian 
Eq. ([3]), and moving within a closed box extending from x = to x = L. The 
particles execute usual Hamiltonian dynamics except when any of the end parti- 
cles hit the walls. Thus when particle I = 1 at the left end (x = 0) hits the wall at 
temperature Tl, it is reflected with a random velocity chosen from the distribution: 

P{v) = e{v) e - miV y(2k B T L ) ; (8) 

where 6(v) is the Heaviside step function. A similar rule is applied at the right end. 

There are two ways of defining a current variable depending on whether one is 
using a discrete or a continuum description. For lattice models, where every particle 
moves around specified lattice points, the discrete definition is appropriate. In a 
fluid system, where the motion of particles is unrestricted, one has to use the 
continuum definition. For the ID hard particle gas, the ordering of particles is 
maintained, and in fact both definitions have been used in simulations to calculate 
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the steady state current. We will show here explicitly that they are equivalent. Let 
us first discuss the discrete definition of heat current. 

For the Langevin and Nose-Hoover baths, we note that the equation of motion 
has the form pi = fi + Si,i/l + ^i,n/r where fi and Jr are forces from the bath. 
The instantaneous rate at which work is done by the left and right reservoirs on 
the system are respectively given by: 

jl,L = fLVl 
and 3n,r = Irvn , 

and these give the instantaneous energy currents from the reservoirs into the sys- 
tem. To define the local energy current inside the wire we first define the local 
energy density associated with the Z th particle (or energy at the lattice site /) as 
follows: 

ei = ^ + V{x l ) + \u{x l -x 2 ) , 
lm\ 2 

ei = -^- + V(xi) + hu(x l . 1 -x l ) + U(xi-xi +1 )], for / = 2,3..JV-1 
2mi 2 

e N = + V(x N ) + ~U{x N -i - x N ) . (9) 

Taking a time derivative of these equations, and after some straightforward ma- 
nipulations, we get the discrete continuity equations given by: 

ei = -32,x + ji,L 

ei = -j/+i,/ + jij-i for I = 2, 3..JV - 1 
cat = 3n,r + 3N,N-l , (10) 

with j tjl _i = 2(^-1 + vi)f 1,1-1 (11) 
and where f ht+l = -fi+ij = -d Xl U(xi - x l+1 ) 

is the force that the (I + l) th particle exerts on the I th particle and v\ = x\. From 
the above equations one can identify ji,i—i to be the energy current from site I — 1 
to I. The steady state average of this current can be written in a slightly different 
form which has a clearer physical meaning. We will denote steady state average 
of any physical quantity A by {A). Using the fact that (dU(xi_i — xi)/dt) = it 
follows that (uj-i/y-i) = (vifij-i) and hence: 

(jl,l-i) = {^{vi + vi-i)fi t i-x) = (vrfij-x) , (12) 

and this has the simple interpretation as the average rate at which the (I — l) th 
particle does work on the I th particle. In the steady state, from Eq. (fTU|) . we get 
the equality of current flowing between any neighbouring pair of particles: 

J = (ji,l) = (32,1} = (33,2) = -(3n,n-i) = -Un,r) , (13) 



where we have used the notation J for the steady state energy current per bond. 
In simulations one can use the above definition, which involves no approximations, 
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and a good check of convergence to steady state is to verify the above equality on 
all bonds. In the case where interaction is in the form of hard particle collisions we 
can write the expression for steady state current in a different form. Replacing the 
steady state average by a time average we get: 

0m-i> = (vifl,l-i) = lim - ( = lim - VAif M _i , (14) 

where t c denotes time instances at which particles I and (I — 1) collide and Aif/ j_i 
is the change in energy of the / th particle as a result of the collision. 

Next we discuss the continuum definition of current which is more appropriate for 
fluids but is of general validity. We will discuss it for the case of Maxwell boundary 
conditions with the system confined in a box of length L. Let us define the local 
energy density at position x and at time t as: 

JV 

e(x,t) = ^2e l 5[x-xi(t)] , (15) 
1=1 

where e/ is as defined in Eq. 0. Taking a time derivative we get the required 
continuum continuity equation in the form (for < x < L ): 

de(x,t) dj(x,t) . . . . . 

— dt dx — =J1 > L °( X > + 3n,rS{x - L) (16) 

where j(x, t) = j K (x, t) + ji(x, t) 

N 

with j K (x,t) = ^eiitfviitfdlx -»/(*)] 
1=1 

N-l 

and ji(x,t) = ^ 0[x-xi(t)] + j 2 ,i0[x - xi(t)] - j N:N -i9[x - x N (t)} . 

1=2 

Here ji t i-i, ji t L, 3n,r are as defined earlier in the discrete case, and we have 
written the current as a sum of two parts, jx and jj, whose physical meaning we 
now discuss. To see this, consider a particle configuration with xx,x 2 , ■■■Xu, < x < 
Xk+i,Xk+2, -XN- Then we get 

3l(x,t) = jk+i,k 

which is thus simply the rate at which the particles on the left of x do work on 
the particles on the right. Hence we can interpret ji(x,t) as the contribution to 
the current density coming from interparticle interactions. The other part jx(x,t) 
arises from the physical flow of particles carrying energy across the point x. Note, 
however, that even in the absence of any net convection particle flow, both jx and 
ji can contribute to the energy flow. In fact for point particles interacting purely 
by hard elastic collisions jk+i k 1S zero whenever the k th and (k + l) th particles are 
on the two sides of the point x and hence ji is exactly zero. The only contribution 
to the energy current then comes from the part jx and we thus for the steady- 
state current we obtain 

N 3 

i=i 
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In simulations either of expressions Eq. (114D or Eq. (jlTf) can be used to evaluate the 
steady state current and will give identical results. For hard particle simulations 
one often uses a simulation which updates between collisions and in this case it is 
more efficient to evaluate the current using Eq. (I14p . We now show that, in the 
nonequilibrium steady state, the average current from the discrete and continuum 
definitions are the same, i. e. , 

= (j(x,t)) = J . (18) 

Note that the steady state current is independent of I or x. To show this we first 
define the total current as: 

J(t)= dxj(x,t)= ^ em - y~] Xi(ji+i t i-ji t i-i)-Xij2,l + XNjN,N-l 
^° l=l,N l=2,N-l 

= ^2 em + ^2 ( xi+i ~ x i)h+i,i ■ ( i9 ) 

l=l,N l=l,N-l 



Taking the steady-state average of the above equation and using the fact that 
(em) = -(eixi) = - ji,i-i)xi), where j 1>0 = ji,L,3N+i,N = -Jn,r we get: 



(J) = -(xijlL + X N j NR ) . 



Since the Maxwell baths are located at x = and x = L, the above then gives 
(J) = L(j(x,t)) = —L(Jnr) and hence from Eq. (fT3j) . we get (j(x,t)) = J = 
(ju—i), which proves the equivalence of the two definitions. 

The extensions of the current definitions, both the discrete and continuum ver- 
sions, to higher dimensions is straightforward. Here, for reference, we outline the 
derivation for the continuum case since it is not easy to see a discussion of this in 
the literature. Consider a system in d-dimensions with Hamiltonian given by: 



H 



E 



(20) 



where x; 



(xj,x 



2 



xf) and pi = (pj ,pf , ...pf) are the vectors denoting the po- 

\xi — x n |. The particles are 
As before we define the local 



sition and momentum of the / th particle and r\ n = 
assumed to be inside a hyper cubic box of volume L d . 
energy density as: 



e(x, t) = ^(x — x/)e/ where 



Taking a derivative with respect to time (and suppressing the source terms arising 
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from the baths) gives: 



de(x,t) 

at 



o=l 



d 



dx al 



+ ^(5(x-X;)q 



(21) 



where 



j£(x,t) = £<5(x-x z ) Q 



"1 



and j?(x,t)=-££0(* a -xf) 

l nj4 



Yl 5{x y - xi) j l>n 



(22) 



where j hn = ^(< + <) /£« » 



and /" n = —dU(r^ n )/dxf is the force, in the a th direction, on the I th particle due 
to the n th particle. We have defined ji n as the current, from particle n to particle 
/, analogously to the discrete ID current. The part jf gives the energy flow as 
a result of physical motion of particles across x a . The part jf also has a simple 
physical interpretation, as in the ID case. First note that we need to sum over 
only those n for which < x a . Then the formula basically gives us the net rate, 
at which work is done, by particles on the left of x a , on the particles to the right. 
This is thus the rate at which energy flows from left to right. By integrating the 
current density over the full volume of the system, we get the total current: 



*n*)=£^f +-££(; 



(23) 



Thus we get an expression similar to that in ID given by Eq. (|19|) . In simulations 
making nonequilibrium measurements, any of the various definitions for current 
can be used to find the steady state current. However, it is not clear whether other 
quantities, such as correlation functions obtained from the discrete and continuum 
definitions, will be the same. 

The local temperature can also be defined using either a discrete approach (giv- 
ing Ti) or a continuum approach (giving T(x,t)). In the steady state these are 
respectively given by (in the ID case): 



k B T(x) 



El 

mi 



*0> 



(24) 



Again it is not obvious that these two definitions will always agree. Lattice simu- 
lations usually use the discrete definition while hard particle simulations use the 
continuum definition. 
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The precise definition of thermal conductivity would be: 

k = lim lim — — , (25) 
L^co AT^o AT y ' 

where AT = Tl — Tr. In general k would depend on temperature T. The finiteness 
of k(T), along with Fourier's law, implies that even for arbitrary fixed values of 
Tl,Tr the current J would scale as ~ (or N^ 1 ). What is of real interest is 
this scaling property of J with system size. We will typically be interested in the 
large N behaviour of the conductivity defined as: 

JN . . 

Kn = AT ' ^ ' 

which will usually be denoted by k. For large N, systems with normal diffusive 
transport give a finite k while anomalous transport refers to the scaling 

k ~ N a tt^O, (27) 

and the value of the heat conduction exponent a is one of the main objects of 
interest. For the current, this implies J ~ iV" -1 . 

The only examples where the steady state current can be analytically evaluated, 
and exact results are available for the exponent a, corresponds to harmonic lattices, 
using very specific methods that will be discussed in sec. ([3]). 

Coupling to baths and contact resistance: In the various models of heat 
baths that we have discussed, the efficiency with which heat exchange takes place 
between reservoirs and system depends on the strength of coupling constants. For 
example, for the Langevin and Nose-Hoover baths, the parameters 7 and 9 re- 
spectively determine the strength of coupling ( for the Maxwell bath one could 
introduce a parameter which gives the probability that after a collision the parti- 
cle's speed changes and this can be used to tune the coupling between system and 
reservoir). From simulations it is found that typically there is an optimum value of 
the coupling parameter for which energy exchange takes place most efficiently, and 
at this value one gets the maximum current for given system and fixed bath tem- 
peratures. For too high or too small values of the coupling strength the current is 
small. The coupling to bath can be thought of as giving rise to a contact resistance. 
An effect of this resistance is to give rise to boundary jumps in the temperature 
profile measured in simulations. One expects that these jumps will be present as 
long as the contact resistance is comparable to the systems resistance. We will later 
see that in order to be sure that one is measuring the true resistance of the system, 
it is necessary to be in parameter regimes where the contact resistances can be 
neglected. 



2.2. Green- Kubo Formula 

The Green-Kubo formula provides a relation between transport coefficients, such 
as the thermal conductivity n or the electrical conductivity a, and equilibrium time 
correlation functions of the corresponding current. For the thermal conductivity in 
a classical ID system, the Green-Kubo formula gives: 



lim lim - ( dt{J(0)J{t)) , (28) 
■— >oo L^oo L Jq 
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where J is the total current as defined in Eq. (|19h and (...) denotes an average over 
initial conditions chosen from either a micro-canonical ensemble or a canonical one 
at temperature T. Two important points to be remembered with regard to use of 
the above Green-Kubo formula are the following: 

(i) It is often necessary to subtract a convective part from the current definition 
or, alternatively, in the microcanonical case one can work with initial conditions 
chosen such that the centre of mass velocity is zero (see also discussions in 0, S] ) . 
To understand this point, let us consider the case with V(x) = 0. Then for a closed 
system that is not in contact with reservoirs, we expect the time average of the 
total current to vanish. But this is true only if we are in the centre of mass frame. 
If the centre of mass is moving with velocity v then the average velocity of any 
particle (vj) = v. Transforming to the moving frame let us write V; = v^ + v. Then 
the average total current in the rest frame is given by (in d-dimensions): 



M 2 , V-/ /v 



(E + PV) v a , 



+E 



(29) 



where M = Yl m i an d E is the average total energy of the system as measured in 
the rest frame and V = L d . In deriving the above result we have used the standard 
expression for equilibrium stress-tensor given by: 



Va av = ^imv'fv'n + ^ E< ( x i ~ <) fin ) , 



(30) 



and assumed an isotropic medium. Thus, in general, to get the true energy current 
in an arbitrary equilibrium ensemble one should use the expression: 



J c a = J*-(E + PV)v c 



(31) 



The corresponding form in ID should be used to replace J in Eq. (|28p . 

(ii) The second point to note is that in Eq. (|28p the order of limits L — > oo and 
then t — ► oo has to be strictly maintained. In fact for a system of particles inside 
a finite box of length L it can be shown exactly that: 



dt{Jc(P)Jc(t)) = . (32) 



To prove this, let us consider a microcanonical ensemble (with {J) = 0, so that 
J c = J), in which case from Eq. (|19p we get: 



J(1 » - it 



N 



i=i 



(33) 



Multiplying both sides of the above equation by J7"(0), integrating over t and noting 
that both the boundary terms on the right hand side vanish, we get the required 
result in Eq. (|32p . With the correct order of limits in Eq. (|28p , one can calculate the 
correlation functions with arbitrary boundary conditions and apply the formula to 
obtain the response of an open system with reservoirs at the ends. 
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Derivation of the Green-Kubo formula for thermal conductivity: There 
have been a number of derivations of this formula by various authors including 
Green, Kubo, Mori, McLennan, Kadanoff and Martin, Luttinger, and Visscher [4, 
0, 0, 0, @, S E3, 11, 12] ■ None of these derivations are rigorous and all require certain 



assumptions. Luttinger's derivation was an attempt at a mechanical derivation and 
involves introducing a fictitious 'gravitational field' which couples to the energy 
density and drives an energy current. However one now has to relate the response 
of the system to the field and its response to imposed temperature gradients. This 
requires additional inputs such as use of the Einstein relation relating diffusion 
coefficient (or thermal conductivity) to the response to the gravitational field. We 
believe that this derivation, as is also the case for most other derivations , implicitly 
assumes local thermal equilibrium. Although these derivations are not rigorous, 
they are quite plausible, and it is likely that the assumptions made are satisfied in 
a large number of cases of practical interest. Thus the wide use of the Green-Kubo 
formula in calculating thermal conductivity and transport properties of different 
systems is possibly justified in many situations. 

Here we give an outline of a non-mechanical derivation of the Green-Kubo for- 
mula, one in which the assumptions can be somewhat clearly stated and their 
physical basis understood. The assumptions we will make here are: 

(a) The nonequilibrium state with energy flowing in the system can be described 
by coarse-grained variables and the condition of local thermal equilibrium is satis- 
fied. 

(b) There is no particle flow, and energy current is equal to heat current. 
The energy current satisfies Fourier's law which we write in the form J(x, t) = 
—Ddu(x,t)/dx where D = k/c, c is the specific heat capacity, and u(x,t) = 
e(x,t), J(x,t) = j(x,t) are macroscopic variables obtained by a coarse-graining 
(indicated by bars) of the microscopic fields. 

(c) Finally, we assume that regression of equilibrium energy fluctuations occurs 
in the same way as nonequilibrium flow of energy. 

We consider a macroscopic system of size L. Fluctuations in energy density in 
equilibrium can be described by the correlation function S(x, t) = (e(x, t)e(0, 0)) — 
(e(x,t)) (e(0, 0)). Assumption (c) above means that the decay of these fluctuations 
is determined by the heat diffusion equation and given by: 

dt dx 2 t>u ' 

where we have also assumed temperature fluctuations to be small enough so that 
the temperature dependence of D can be neglected. From time reversal invariance 
we have S(x,t) = S(x, —t). Using this and the above equation we get: 

„ wt = 2Dk 2 S(k,t = 0) 
D 2 k 4 + to 2 ' 



S(k,u) = / dtS(k,t)e 



where S(k,t) = dxS(x,t)e lkx . Now, from equilibrium statistical mechanics we 
have S(k = 0, t = 0) = cksT 2 and using this in the above equation we obtain: 



1 .. .. u? 



K = cD = ^SSf*^ • (34) 



One can relate the energy correlator to the current correlator using the continuity 
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equation de(x,t)/dt + dj(x,t)/dx = and this gives: 

Q{k,oj)~j{k' = (27r) 2 ^S(k,io) 5(k + k')S(u; + u') , 

where j(k,u>) = J^° oo dxdtj(x,t)e l( - kx ^ ujt \ Integrating the above equation over k! 
gives: 

oo poo 2 

dx / dt <j(x,t)j(0,0)) e< kx ~^ = —S(k,co) . (35) 

-oo J —oo ™ 



From Eq. (|34l) and Eq. (|35j) we then get: 

-i poo poo 

k = lim lim — — dx dt <jOM)j(0,0)) e^~^ . (36) 
^Ok^o2k B T z J^ce 

Finally, using time-reversal and translational invariance and interpreting the 
lo —* 0, k —* limits in the alternative way ( as r — > oo, L — > oo) we recover the 
Green-Kubo formula in Eq. (|28p . 

Limitations on use of the Green-Kubo formula: There are several situa- 
tions where the Green-Kubo formula in Eq. (|28p is not applicable. For example, 
for the small structures that are studied in mesoscopic physics, the thermodynamic 
limit is meaningless, and one is interested in the conductance of a specific finite sys- 
tem. Secondly, in many low dimensional systems, heat transport is anomalous and 
the thermal conductivity diverges. In such cases it is impossible to take the limits 
as in Eq. (|28|) : one is there interested in the thermal conductance as a function of 
L instead of an L-independent thermal conductivity. The usual procedure that has 
been followed in the heat conduction literature is to put a cut-off at t c ~ L, in the 
upper limit in the Green-Kubo integral 0]. The argument is that for a finite system 
connected to reservoirs, sound waves traveling to the boundaries at a finite speed, 
say v, lead to a decay of correlations in a time ~ L/v. However there is no rigorous 
justification of this assumption. A related case is that of integrable systems, where 
the infinite time limit of the correlation function in Eq. (I28p is non-zero. 

Another way of using the Green-Kubo formula for finite systems is to include 
the infinite reservoirs also while applying the formula and this was done, for ex- 



ample, by Allen and Ford [13] for heat transport and by Fisher and Lee LJ] for 
electron transport. Both these cases are for non-interacting systems and the final 
expression for conductance (which is more relevant than conductivity in such sys- 
tems) is basically what one also obtains from the Landauer formalism [15| |. or the 
nonequilibrium Green's function approach [see sec. (|2.3p ]. 

It has been shown that Green-Kubo like expressions for the linear response heat 
current for finite open systems can be derived rigorously by using the steady state 



fluctuation theorem 116, 17, 18, 19, 20, 21, 221. This has been done for lattice mod 



els coupled to stochastic Markovian baths and the expression for linear response 
conductance of a one dimensional chain is given as: 

T -I pOO 

lim-^- = — =/ (ji(0)m), (37) 



AT^o AT k B T 2 

where ji is the discrete current defined in sec. (12.1(1 . Some of the important dif- 
ferences of this formula with the usual Green-Kubo formula are worth keeping in 
mind: (i) the dynamics of the system here is non-Hamiltonian since they are for a 
system coupled to reservoirs, (ii) one does not need to take the limit iV — > oo first, 
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the formula being valid for a finite system, (iii) the discrete bond current appears 
here, unlike the continuum one in the usual Green-Kubo formula. Recently, an ex- 
act linear response result similar to Eq. (137p . for the conductance of a finite open 
system has been derived using a different approach [23]. This has been done for 
quite general classical Hamiltonian systems and for a number of implementations 
of heat baths. 

It appears that the linear response formula given by Eq. (|37j) is the correct 
one to use to evaluate the conductance in systems where there is a problem with 
the usual Green-Kubo formula, e.g. in finite systems or low dimensional systems 
showing anomalous transport (because of slow decay of {J{Q)J{t)). We note here 
the important point that the current-current correlation can have very different 
scaling properties, for a purely Hamiltonian dynamics, as compared to a heat bath 



dynamics. This has been seen in simulations by Deutsche and Narayan [2J] for 
the random collision model [defined in sec. ([4.2. ip ]. Of course this makes things 
somewhat complicated since the usefulness of the Green-Kubo formalism arises 
from the fact that it allows a calculation of transport properties from equilibrium 
properties of the system, and without any reference to heat baths, etc. In fact, as 
we will see in sec. (|4.ip . all the analytic results for heat conduction in interacting 
one dimensional systems rely on Eq. (|28p and involve a calculation of the current- 
current correlator for a closed system. Some of the simulation results discussed later 
suggest that, for interacting (nonlinearly) systems, in the limit of large system size 
the heat current is independent of details of the heat baths. This means that, in 
the linear response regime and the limit of large system sizes, a description which 
does not take into account bath properties may still be possible. 



2.3. Nonequilibrium Green's function method 

The nonequilibrium Green's function method (NEGF) is a method, first invented 
in the context of electron transport, to calculate steady state properties of a finite 
system connected to reservoirs which are themselves modeled by noninteracting 



Hamilitonians with infinite degrees of freedom [46|, 1471 . |48| |. Using the Keldysh 
formalism, one can obtain formal expressions for the current and other observ- 
ables such as electron density, in models of electrons such as those described by 
tight-bindin g ty pe Hamiltonians. Recently this approach has been applied both to 
phonon 50, 51] and photon 55|, 5(| transport. 



The main idea in the formalism is as follows. One starts with an initial density 
matrix describing the decoupled system, and two infinite reservoirs which are in 
thermal equilibrium, at different temperatures and chemical potentials. The system 
and reservoirs are then coupled together and the density matrix is evolved with the 
full Hamiltonian for an infinite time so that one eventually reaches a nonequilibrium 
steady state. Various quantitites of interest such as currents and local densities, etc 
can be obtained using the steady state density matrix and can be written in terms of 
the so called Keldysh Green's functions. An alternative and equivalent formulation 
is the Langevin approach where, instead of dealing with the density matrix, or 
in the classical case with the phase space density, one works with equations of 
motion of the phase space variables of the full Hamiltonian (system plus reservoirs) . 
Again the reservoirs, which are initially in thermal equilibrium, are coupled to the 
system in the remote past. It is then possible to integrate out the reservoir degrees 
of freedom, and these give rise to generalized Langevin terms in the equation of 
motion. For non-interacting systems, one can show that it is possible to recover all 
the results of NEGF exactly, both for electrons and phonons. Here we will discuss 
this approach for the case of phonons, and describe the main results that have been 
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obtained [sec. (13, 2j) ]. 

For non- interacting systems, the formal expressions for current obtained from 
the NEGF approach is in terms of transmission coefficients of the heat carriers 
(electrons, phonons or photons) across the system, with appropriate weight fac- 
tors corresponding to the population of modes in the reservoirs. These expressions 



are basically what one also obtains from the Landauer formalism 151 ]. We note 
that in the Landauer approach one simply thinks of transport as a transmission 
problem and the current across the system is obtained directly using this picture. 
In the simplest set-up one thinks of one-dimensional reservoirs (leads) filled with 
non-interacting electrons at different chemical potentials. On connecting the sys- 
tem in between the reservoirs, electrons are transmitted through the system from 
one reservoir to the other. The net current in the system is then the sum of the 
currents from left-moving and right-moving electron states from the two reservoirs 
respectively. 



3. Heat conduction in harmonic lattices 

The harmonic crystal is a good starting point for understanding heat transport 
in solids. Indeed in the equilibrium case we know that studying the harmonic 
crystal already gives a good understanding of, for example, the specific heat of 
an (electrically) insulating solid. In the nonequilibrium case, the problem of heat 
conduction in a classical one dimensional harmonic crystal was studied for the 
first time by Rieder, Lebowitz and Lieb (RLL) [2^]. They considered the case of 
stochastic Markovian baths and were able to obtain the steady state exactly. The 
main results of this paper were: (i) the temperature in the bulk of the system was 
constant and equal to the mean of the two bath temperatures, (ii) the heat cur- 
rent approaches a constant value for large system sizes and an exact expression for 
this was obtained. These results can be understood physically when one realizes 
that in the ordered crystal, heat is carried by freely propagating phonons. RLL 
considered the case where only nearest neighbour interparticle interactions were 
present. Nakazawa (NK) [2^] extended these results to the case with a constant 
onsite harmonic potential at all sites and also to higher dimensions. The approach 
followed in both the RLL and NK papers was to obtain the exact nonequilibrium 
stationary state measure which, for this quadratic problem, is a Gaussian distri- 
bution. A complete solution for the correlation matrix was obtained and from this 
one could obtain both the steady state temperature profile and the heat current. 

In sec. (]3.ip we will briefly describe the RLL formalism. In sec. (]3.2[) we will 
describe a different and more powerful formalism. This is the Langevin equation 
and nonequilibrium Green's function method and various applications of this will 
be discussed in sees. (|3.3l3,4l3,5p . 



3.1. The Rieder- Lebowitz- Lieb method 

Let us consider a classical harmonic system of N particles with displacements 
about the equilibrium positions described by the vector X = {x\, X2, ...xn} t , where 
T denotes transpose. The particles can have arbitrary masses and we define a 
diagonal matrix M whose diagonal elements Mu = mi, for I = 1,2...N, give the 
masses of the particles. The momenta of the particles are given by the vector 
P = (MX) = {pi,P2, ■■■Pn} T • We consider the following harmonic Hamiltonian for 
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the system: 



H = -P T M~ 1 P + -X T $X , (38) 



where is the force matrix. Let us consider the general case where the I th particle is 
coupled to a white noise heat reservoir at temperature T B with a coupling constant 
7;. The equations of motion are given by: 

Pi 

xi = — 
mi 

Pl = -Y,®lnX n -—Pl + Vl for | = 1,2...JV (39) 

n 1 

with the noise terms satisfying the usual fluctuation dissipation relations 

(m(t)Vn(t')} = 2 ll k B T l B 5 l , n 5{t - t') . (40) 

Defining new variables q = {qi,q2---Q2N} T = {^l, X2---Xn,Pi,P2---Pn} T we can 
rewrite Eqs. (|39I40|) in the form: 

q = —A q + r] 

mi 1 T {t')) = b5(t-t i ) (4i) 

where the 2N dimensional vector rj T = (0, 0, ..0, T]i,t]2---VN) and the 2N x 2N 
matrices A, D are given by: 



,0 -M- x \ ~ _ (0 

1 $M-!f ) V -\0E 



and f and E are N x N matrices with elements T[ n = jiSi n , E( n = 2kBT l B ^/i5i n 
respectively. In the steady state, time averages of total time derivatives vanish, 
hence we have (d(qq T ) / dt) = 0. From Eq. (I41D we get 

(±(qq T )) = ((-Aq + V )q T + q(-q T A T + rf)) 

= -A.B - B.A T + (rjq T + qr] T ) = (42) 



where B is the correlation matrix (qq T ). To find the average of the term involving 
noise we first write the formal solution of Eq. [3TJ 



q(t) = G{t - t )q(t ) + [ dt'G{t - t')r]{t') where G(t) = e~ At . (43) 

Jtn 
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Setting to — ► — °° and assuming G(oo) = (so that a unique steady state exists), 
we get q(t) = j 1 ^ dt'G(t - t')r/(t') and hence 



(q V T )= \ dt'Git-t'Xvit'^it)) 
= / dt'G(t - t')D5{t - t') = -D (44) 

J-oo ^ 



where we have used the noise correlations given by Eq. (|40p . and the fact that 
G(0) = I, a unit matrix. Using Eq. ([44]) in Eq. (|42|) . we finally get 



+ = L> (45) 

The solution of this equation gives the steady state correlation matrix B which 
completely determines the steady state since we are dealing with a Gaussian pro- 
cess. In fact the steady state is given by the Gaussian distribution 

P({ qi }) = (2tt)- n Det[B\- ll2 e~^ qT ^ lq . (46) 



Some of the components of the matrix equation Eq. (|45l) have simple physical 
interpretations. To see this we first write B in the form 



B 



B X p 



where B x , B p and B xp are NxN matrices with elements (B x )i n = (xix n ), (B p )i n = 
(pip n ) and (B xp )i n = (xip n ). From Eq. (f4"5j) we then get the set of equations: 

M- X B T £V + B xp M~ l = (47) 
B X $ T - AT 1 B p + BxpM^f = (48) 
M _1 f B p + BpfM- 1 + $B xp + Blp$ = E (49) 

From Eq. H7]we get the identity {xip n /m n ) = —(x n pi/mi) . Thus (xipi) = 0. The 
diagonal terms in Eq. 08] give 



V(x,(*i„x„)) - — (pf) + {xipi)— = 

— * m j m i 



mi ffi\ 

* (xffi = (&) = k B l\ 
oxi mi 

where we have defined the local temperature k B T[ = (pf/mi). This equation has 
the form of the 'equipartition' theorem of equilibrium physics. It is in fact valid 
quite generally for any Hamiltonian system at all bulk points and simply follows 
from the fact that ((d/dt)xipi) = 0. Finally let us look at the diagonal elements of 
Eq. 09] This gives the equation 

yVC^Zn)— ) = —k B (Tf - Ti) (50) 
^— ' mi mi 
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which again has a simple interpretation. The right hand can be seen to be equal 
to ((—"tipi/mi + rji)pi/mi) which is simply the work done on the Zth particle by 
the heat bath attached to it (and is thus the heat input at this site). On the left 
hand side {{—<&i n Xn)pi / mi) is the rate at which the n th particle does work on the 
I th particle and is therefore the energy current from site n to site I. The left side is 
thus is just the sum of the outgoing energy currents from the Ith site to all the sites 
connected to it by Thus we can interpret Eq. (I50j) as an energy-conservation 
equation. The energy current between two sites is given by 



Our main interests are in computing the temperature profile and the energy current 
in the steady state. This requires solution of Eq. ([13]) and this is quite difficult and 
has been achieved only in a few special cases. For the one-dimensional ordered 
harmonic lattice, RLL were able to solve the equation exactly and obtain both 
the temperature profile and the current. The extension of their solution to higher 
dimensional lattices is straightforward and was done by Nakazawa. More recently 
Bonetto et al [(33] have used this approach to solve the case with self-consistent 
reservoirs attached at all the bulk sites of a ordered harmonic lattice. 

A numerical solution of Eq. (|45p requires inversion of a N(2N + 1) x N(2N + 1) 
matrix which restricts one to rather small system sizes. The RLL approach is 
somewhat restrictive since it is not easily generalizable to other kinds of heat 
baths or to the quantum case. Besides, except for the ordered lattice, it is difficult 
to obtain useful analytic results from this approach. In the next section we discuss a 
different approach which is both analytically more tractable, as well as numerically 
more powerful. 

3.2. Langevin equations and Green's function (LEGF) formalism 

This approach involves a direct solution of generalized Langevin equations. Using 
this solution one can evaluate various quantities of interest such as steady state 
current and temperature profiles. Compact formal expressions for various quanti- 
ties of interest can be obtained and, as pointed out earlier, it turns out that these 
are identical to results obtained by the nonequilibrium Green's function (NEGF) 
method described in sec. (|2.3p . The method can be developed for quantum me- 
chanical systems, in which case we deal with quantum Langevin equations (QLE), 
and we will see that the classical results follow in the high temperature limit. In all 
applications we will restrict ourselves to this approach which we will henceforth re- 
fer to as the Langevin equations and Green's function (LEGF) method. As we will 
see, for the ordered case to be discussed in sec. (13.31) . one can recover the standard 
classical results as well as extend them to the quantum domain using the LEGF 
method. For the disordered case too [sec. fl3.4|> ] . one can make significant progress. 
Another important model that has been well studied in the context of harmonic 
systems is the case where self-consistent heat reservoirs are attached at all sites 
of the lattice. In sec. (13.5j) we will review results for this case obtained also from 
the LEGF approach. All examples in this section deal with the case where particle 
displacements are taken to be scalars. The generalization to vector displacements 
is straightforward. 

The LEGF formalism has been developed in the papers by (4ol. l4ll . I42I . I431 . 13] 
and relies on the approach first proposed by Ford, Kac and Mazur [45j of modeling 
a heat bath by an infinite set of oscillators in thermal equilibrium. Here we will 
outline the basic steps as given in (iij. One again considers a harmonic system 




(51) 
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which is coupled to reservoirs which are of a more general form than the white 
noise reservoirs studied in the last section. The reservoirs are now themselves taken 
to be a collection of harmonic oscillators, whose number will be eventually taken 
to be infinite. As we will see, this is equivalent to considering generalized Langevin 
equations where the noise is still Gaussian but in general can be correlated. We 
will present the discussion for the quantum case and obtain the classical result as 
a limiting case. 

We consider here the case of two reservoirs, labeled as L (for left) and R (right) 
, which are at two different temperatures. It is easy to generalize to the case where 
there are more than two reservoirs. For the system let X = {x\,X2, ...xn} t now 
be the set of Heisenberg operators corresponding to the displacements (assumed 
to be scalars) of the iV particles, about equilibrium lattice positions. Similarly let 
Xl and X R refer to position operators of the particles in the left and right baths 
respectively. The left reservoir has Nl particles and the right has N R particles. 
Also let P, Pl-, Pr be the corresponding momentum variables satisfying usual com- 
mutation relations with the position operators {i.e. [x;,p n ] = iM/ jn , etc.). The 
Hamiltonian of the entire system and reservoirs is taken to be: 

H = H S + H L + H R + H I L + H I R (52) 
where H s = ^P T M^P + ^X T $X , 

H L = l -PlMl l P L + \x T L ^ L X L , 

H R = l -PlM R l P R + \x T R ® R X R , 
H L = X T V L X L , H R = X T V R X R , 



where M, Ml , M R are real diagonal matrices representing masses of the particles in 
the system, left, and right reservoirs respectively. The quadratic potential energies 
are given by the real symmetric matrices 4*^, & R while Vl and V R denote the 
interaction between the system and the two reservoirs respectively. It is assumed 
that at time t = to, the system and reservoirs are decoupled and the reservoirs are 
in thermal equilibrium at temperatures Tl and T R respectively. 
The Heisenberg equations of motion for the system (for t > to) are: 

MX = —&X — V L X L - V R X R , (53) 

and the equations of motion for the two reservoirs are 

M L X L = -§ L X L - V£X , (54) 
M R X R = -§ R X R - V£X . (55) 



One first solves the reservoir equations by considering them as linear inhomoge- 
neous equations. Thus for the left reservoir the general solution to Eq. (|54p is (for 
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t > t ): 

X L {t) = f£(t - t )M L X L (t ) + g+(t - t )M L X L (t ) 

- fdt' g+(t-t')\?X(t') , (56) 

Jt a 

with f+(t) = U L cos {(l L t)Ul 9(t), g+(t) = U L Slli f Lt) U T L 6(t) , 

where 6(t) is the Heaviside function, and Ul, SIl are the normal mode eigenvector 
and eigenvalue matrices respectively, corresponding to the left reservoir Hamilto- 
nian Hl, and which satisfy the equations: 

Ul&lUl = &l , UlM L U L = I . 

A similar solution is obtained for the right reservoir. Plugging these solutions back 
into the equation of motion for the system, Eq. (|53p . one gets the following effective 
equations of motion for the system: 

MX = -$X + VL + f dt't L (t-t')X(t')+r] R + f dt't R (t-t')X(t'), (57) 

J to J to 

where t L (t) = V L g+(t) V[ , t R (t) = V R g+(t) \% 

and 7] L = -V L [ ft{t ~ t )M L X L (t ) + g+(t - t )M L X L (t ) 

m = -V R \ /£(* - t )M R X R {t ) + g+{t - t )M R X R (t ) 



This equation has the form of a generalized quantum Langevin equation. The 
properties of the noise terms rji and rj R are determined using the condition that, at 
time to, the two isolated reservoirs are described by equilibrium phonon distribution 
functions. At time to, the left reservoir is in equilibrium at temperature Tl and 
the population of the normal modes (of the isolated left reservoir) is given by the 
distribution function /ft(o;,Ti) = l/^ kJ ' kBTL — 1]. One then gets the following 
correlations for the left reservoir noise: 



(VL(t)r)l(t')} = V L U L 



cos Clz(t — t')—^ — coth (- 



"2k B T L ' 



— ism Ox (i — t ) 



2Q, 



mi , 



(58) 



and a similar expression for the right reservoir. The limits of infinite reser- 
voir sizes (Nl,N r — > oo) and to — ► — °° are now taken. One can then solve 
Eq. (|57p by taking Fourier transforms. Let us define the Fourier transforms 
X£) = (1/27T) fr^dt X(t)e^, fj L>R (u) = (1/2*) f^dt r} L;R (t)e^\ g+ R (u) = 
J^oo dt g\ R (t)e lujt . One then gets the following stationary solution to the equations 
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of motion Eq. (f57|) : 

/CO 
du;X(oj)e- iujt , (59) 

with X(u) = G + (cj) [fj L (u) + fj R (u)] , 



where G + (uj) 



1 



[-uj 2 M + $ - E£(w) - E+(u)] 
and E+( W ) = Vtg+( W )V L T , E+H = V^+( W )Vj . 

For the reservoirs one obtains [using Eq. (|56p ] 

-V L X L (u) = fj L (u) + £+X{u;) , 

-VrX r {u) = fj R (u) + t+X{uj) . (60) 



The noise correlations, in the frequency domain, can be obtained from Eq. (|58p 
and we get (for the left reservoir): 

{f, L {u)f L {J))=5{u + u') T l (uj) -[l + MuM] (61) 

where Tl(w) = 7m[Ej(o;)] 

which is a fluctuation-dissipation relation. This also leads to the more commonly 
used correlation: 

i( fj L (u)vl(u;')+fj L (u;')fjl(u;) ) = S(u + u') f L (u>) A ooth (_^_). ( 6 2) 

Similar relations hold for the noise from the right reservoir. The identification 
of G + {uj) as a phonon Green function, with r (uj) as self energy contributions 
coming from the baths, is the main step that enables a comparison of results derived 
by the LEGF approach with those obtained from the NEGF method. This has been 
demonstrated in refn. [44j]. 

Steady state properties: The simplest way to evaluate the steady state current 
is to evaluate the following expectation value for current from left reservoir into 
the system: 

J = -{X T V L X L ). (63) 

This is just the rate at which the left reservoir does work on the wire. Using the 
solution in Eq. (|59l6UI6ip one obtains, after some manipulations: 

1 r°° 

J = — / dcu T{u)hu> [f(u, T L ) - f(u, T R )} . (64) 

where T(u) = 4 Tr[ G + {u) Y L {u))G'[u) V R {u))] , 

and G~(u>) = G + '(uj). This expression for current is of identical form as the NEGF 
expression for electron current (see for example 0, 47, 4p . l49l| 1 and has also been 
derived for phonons using the NEGF approach in ref ns . [50. 1 5 1J| . In fact this expres- 
sion was first proposed by Angelescu et al. and by Rego and Kirczenow [Hj] for 
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a ID channel and they obtained this using the Landauer approach. Their result was 
obtained more systematically later by Blencowe [13]. Note that Eq. (164D above can 
also be written as an integral over only positive frequencies using the fact that the 
integrand is an even function of uj. The factor T(u) is the transmission coefficient 
of phonons at frequency uj through the system, from the left to right reservoir. The 
usual Landauer result for a ID channel precisely corresponds to Rubin's model of 
bath, to be discussed in sec. (I3.4.ip . 

For small temperature differences AT = T^ — Tr « T, where T = (Tj, + Tr)/2, 
i.e. , in the linear response regime the above expression reduces to: 

AT f°° df(uj,T) 
J=—j^T(u)fu,—^. (65) 

(66) 

The classical limit is obtained by taking the high temperature limit huj/kBT — > 0. 
This gives: 



j _ k B AT 
4tt 



/CO 
duj T(u) . (67) 
-co 



One can similarly evaluate various other quantities such as velocity-velocity corre- 
lations and position-velocity correlations. The expressions for these in the classical 
case are respectively: 



K = {XX 



k ' U ' L 1 " du uj G + {u)T l {u)G-{uj) + M« / du uj G + {uj)Tr{uj)G-{uj) , 



C = (XX 1 ) 



lkjJ ' L 1 * duo G + {uj)T L (uj)G-{uj) + !°° du G + {uj)T R {uj)G-{uj) 



The correlation functions K can be used to define the local energy density which 
can in turn be used to define the temperature profile in the non-equilibrium steady 
state of the wire. Also we note that the correlations C give the local heat current 
density. Sometimes it is more convenient to evaluate the total steady state current 
from this expression rather than the one in Eq. (|64p . 

For one-dimensional wires the above results can be shown [32| to lead to ex- 
pressions for current and temperature profiles obtained in earlier studies of heat 
conduction in disordered harmonic chains [HI, H3, E2] ■ 

In our derivation of the LEGF results we have implicitly assumed that a unique 
steady state will be reached. One of the necessary conditions for this is that no 
modes outside the bath spectrum are generated for the combined model of system 
and baths. These modes, when they exist, are localized near the system and any 
initial excitation of the mode is unable to decay. This has been demonstrated and 



discussed in detail in the electronic context [491 ]. 

We note that unlike other approaches such as the Green-Kubo formalism and 
Boltzmann equation approach, the Langevin equation approach explicitly includes 
the reservoirs. The Langevin equation is physically appealing since it gives a nice 
picture of the reservoirs as sources of noise and dissipation. Also just as the Lan- 
dauer formalism and NEGF have been extremely useful in understanding electron 
transport in mesoscopic systems it is likely that a similar description will be useful 
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for the case of heat transport in (electrically) insulating nanotubes, nanowires, etc. 
The LEGF approach has some advantages over NEGF. For example, in the clas- 
sical case, it is easy to write Langevin equations for nonlinear systems and study 
them numerically. Unfortunately, in the quantum case, one does not yet know how 
to achieve this, and understanding steady state transport in interacting quantum 
systems is an important open problem. 



3.3. Ordered harmonic lattices 

As mentioned above, heat conduction in the ordered harmonic chain was first 
studied in the Rieder, Lebowitz and Lieb (RLL) paper and its higher dimensional 
generalization was obtained by Nakazawa. The approach followed in both the RLL 
and Nakazawa papers was to obtain the exact nonequilibrium stationary state 
measure which, for this quadratic problem, is a Gaussian distribution. A complete 
solution for the correlation matrix was obtained and from this one could obtain 
both the steady state temperature profile and the heat current. 



Here we follow refn. 57j to show how the LEGF method, discussed in the previous 
section, can be used to calculate the heat current in ordered harmonic lattices 
connected to Ohmic reservoirs (for a classical system this is white noise Langevin 
dynamics). We will see how exact expressions for the asymptotic current (N — * oo) 
can be obtained from this approach. We also briefly discuss the model in the 
quantum regime and extensions to higher dimensions. 

3.3.1. One dimensional case 

The model considered in [57j is a slightly generalized version of those studied by 
RLL and Nakazawa. An external potential is present at all sites and the pinning 
strength at the boundary sites are taken to be different from those at the bulk 
sites. Thus both the RLL and Nakazawa results can be obtained as limiting cases. 
Also it seems that this model more closely mimics the experimental situation. In 
experiments the boundary sites would be interacting with fixed reservoirs, and the 
coupling to those can be modeled by an effective spring constant that is expected 
to be different from the interparticle spring constant in the bulk. We also note here 
that the constant onsite potential present along the wire relates to experimental 
situations such as that of heat transport in a nanowire attached to a substrate or, 
in the two-dimensional case, a monolayer film on a substrate. Another example 
would be the heat current contribution from the optical modes of a polar crystal. 

Consider N particles of equal masses m connected to each other by harmonic 
springs of equal spring constants k. The particles are also pinned by onsite quadratic 
potentials with strengths k Q at all sites except the boundary sites where the pinning 
strengths are k Q + k' . The Hamiltonian is thus: 



N 2 i N-l 

H = Df^ + + E 5*0=1+1 " x l? + 5 k '( x i + x n) » (68) 



1=1 



where xi denotes the displacement of the I particle from its equilibrium position. 
The particles 1 and N at the two ends are connected to heat baths at temperature 
Tl and Tr respectively, assumed to be modeled by Langevin equations correspond- 
ing to Ohmic baths (S(cj) = ijuj). In the classical case the steady state heat current 
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from left to right reservoir can be obtained from Eq. (167j) and given by [29L I32I]: 

J= k B (T L -T R ) r* ^ 

where T N (u) = 4 T 2 (lj)\G 1n (lu)\ 2 , G{lu) = Z~ x /k 
and Z = [-mL0 2 I + § - ±{uj)}/k , 

where I is a unit matrix, <I> is the force matrix corresponding to the Hamiltonian in 
Eq. (|68p . The N x N matrix £ has mostly zero elements except for En = Ti^n = 
iT(uj) where T(uj) = jui. The matrix Z is tri-diagonal matrix with Z\\ = Z^n = 
(k+k a +k' —mLo 2 —i^fio)/k, all other diagonal elements equal to 2+k a /k— mui 2 /k and 
all off-diagonal elements equal to —1. Then it can be shown easily that |Gijv(w)| = 
l/(k |Ajv|) where A^r is the determinant of the matrix Z. This can be obtained 
exactly. For large N, only phonons within the spectral band of the system can 
transmit, and the integral over u in Eq. (|69p can be converted to one over q to 
give: 

2 7 2 k B (T L -T R ) r du u 2 
k tt J dq \A N \ Z 

with mui 2 = k a + 2k[l — cos (g)]. Now using the result: 

ite. I " T+^s^m = h "" n-arn^ • (71) 

where g\{q) and g2(q) are any two well-behaved functions, one can show that in 
the limit N — > 00, Eq. (|70p gives 



J = jk 2 k B (n-T R) ^—^ 

where A = 2fc(A; - fe') + A:' 2 + + 2fc )^ and n = 2 fc(A; - fc') + ^ . 

m m 

Two different special cases lead to the RLL and Nakazawa results. First in the case 
of fixed ends and without onsite potentials, i.e. k' = k and k a = 0, we recover the 
RLL result M 



jRLL = kk B {T L -T R ) 



v v I 4 



where v = — =- . (73) 
2 2 V v J 7 2 y ' 



27 

In the other case of free ends, i.e. k' = 0, one gets the Nakazawa result [26| : 
N k-/k B (T L -T R ) 



J 1 



XX 4 



where A = — — — ^ . (74) 
2(mfc + 7 2 ) L 2 2V A J k(mk + 7 2 ) v ; 



In the quantum case, in the linear response regime, Eq. (|65l) and similar manip- 
ulations made above for the N — > 00 limit leads to the following final expression 
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^{Tl-Tr) r sin 2 g 2 2 / hw q \ 

J = 731—7^ / dg- FTZm ^ cosech . ( 75 ) 



for current: 



4nk B mT 2 J "A-Ocosg q \2k B T 
where uj 2 = [k a + 2k(l — cos q)]/m . 

While one cannot perform this integral exactly, numerically it is easy to obtain the 
integral for given parameter values. It is interesting to examine the temperature 
dependence of the conductance J /AT. In the classical case this is independent 
of temperature while one finds that at low temperatures the quantum result is 
completely different. For three different cases one finds, in the low temperature 

regime, the following behaviour: 

'T 3 for k' = k,k o = 
J ~ < ~ T for k' = 0, k Q = (76) 
^ ~ 6 k ^[% bT) for k' = 0, k Q ^ , 

where uj d = (k /m) 1 /' 2 . In studies trying to understand experimental work on 
nanosystems [see sec. (jSj)] , the temperature dependence of the conductance is 
usually derived from the Landauer formula, which corresponds to the Rubin model 
of bath. The temperature dependence will then be different from the above results. 

3.3.2. Higher dimensions 

As shown by Nakazawa [26] the problem of heat conduction in ordered har- 
monic lattices in more than one dimension can be reduced to an effectively one- 
dimensional problem. We will briefly give the arguments here and also give the 
quantum generalization. 

Let us consider a d-dimensional hypercubic lattice with lattice sites labeled by the 
vector n = {n a },a = l,2...d, where each n a takes values from 1 to L a . The total 
number of lattice sites is thus N = L\L<i-..Ld- We assume that heat conduction 
takes place in the a = d direction. Periodic boundary conditions are imposed 
in the remaining d — 1 transverse directions. The Hamiltonian is described by a 
scalar displacement X n and, as in the ID case, we consider nearest neighbour 
harmonic interactions with a spring constant k and harmonic onsite pinning at 
all sites with spring constant k Q . All boundary particles at = 1 and = 
are additionally pinned by harmonic springs with stiffness k' and follow Langevin 
dynamics corresponding to baths at temperatures Tl and Tr respectively. 

Let us write n = (n t ,n^) where n t = (rti, ^■••^d-i)- Also let q = (ffi, 92---9d!-i) 
with q a = 2ns/ L a where s goes from 1 to L a . Then defining variables 

^nd(q) = 1/21/2 7T7T y^^n t ,n d e' q,nt , (77) 
L x L 2 ■■■Li d _ x nt 

one finds that, for each fixed q, X Hd (q) (rid = 1 ; 2...Ld) satisfy Langevin equations 
corresponding to the ID Hamiltonian in Eq. (|68p with the onsite spring constant 
k Q replaced by 



A(q) = K + 2 k [ d - 1 - Y, cos (q a ) ] . (78) 

a=l,d-l 

For Ld — > oo, the heat current J(q) for each mode with given q is then simply 
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given by Eq. (172p with k Q replaced by A q . In the quantum mechanical case we use 
Eq. (j75j) . The heat current per bond is then given by: 

J= T r 1 j E J (q) 

L 1 L 2 ....L ( i-\ ^ 

Note that the result holds for finite lengths in the transverse direction. For infinite 
transverse lengths we get J = / ... J Q 27r <iqJ(q)/(27r) a!-1 . 

3.4. Disordered harmonic lattices 

From the previous section we see that the heat current in an ordered harmonic 
lattice is independent of system size (for large systems) and hence transport is 
ballistic. Of course this is expected since there is no mechanism for scattering 
of the heat carriers, namely the phonons. Two ways of introducing scattering of 
phonons are by introducing disorder in the system, or by including anharmonicity 
which would cause phonon-phonon interactions. In this section we consider the 
effect of disorder on heat conduction in a harmonic system. 

Disorder can be introduced in various ways, for example by making the masses of 
the particles random as would be the case in a isotopically disordered solid, or by 
making the spring constants random. Here we will discuss the case of mass-disorder 
only since the most important features do not seem to vary much with the type of 
disorder one is considering. Specifically, we will consider harmonic systems where 
the mass of each particle is an independent random variable chosen from some 
fixed distribution. 

It can be expected that heat conduction in disordered harmonic systems will 
be strongly affected by the physics of Anderson localization. In fact the problem 
of finding the normal modes of the harmonic lattice can be directly mapped to 
that of finding the eigenstates of an electron in a disordered potential (in a tight- 
binding model, for example) and so we expect the same kind of physics as in 
electron localization. In the electron case the effect of localization is strongest in one 
dimensions where it can be proved rigorously that all eigenstates are exponentially 
localized, hence the current decays exponentially with system size and the system is 
an insulator. This is believed to be true in two dimensions also. In the phonon case 
the picture is much the same except that, in the absence of an external potential, 
the translational invariance of the problem leads to the fact that low frequency 
modes are not localized and are effective in transporting energy. Another important 
difference between the electron and phonon problems is that electron transport 
is dominated by electrons near the Fermi level while in the case of phonons, all 
frequencies participate in transport. These two differences lead to the fact that the 
disordered harmonic crystal in one and two dimensions is not a heat insulator, 
unlike its electronic counterpart. Here we will present results using the LEGF 
approach to determine the system size dependence of the current in one dimensional 
mass-disordered chains. We note that basically this same approach (NEGF) is also 
popular in the electron case and is widely used in mesoscopic physics. Also earlier 
treatments by, for example, Rubin and Greer [13] and by Casher and Lebowitz [2^] 
of the disordered harmonic chain, can be viewed as special examples of the LEGF 
approach. We will also discuss results of simulations for the two-dimensional case. 

Our main conclusions here will be that Fourier's law is not valid in a disordered 
harmonic crystal in one and two dimensions, the current decays as a power law 
with system size and the exponent a is sensitive to boundary conditions (BC) and 
spectral properties of the heat baths. 



(79) 
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3.4-.1. One dimensional disordered lattice 
We first briefly review earlier work on this problem [2] . The thermal conductivity 



of disordered harmonic lattices was first investigated by Allen and Ford [13| who, 
using the Kubo formalism, obtained an exact expression for the thermal conduc- 
tivity of a finite chain attached to infinite reservoirs. From this expression they 
concluded, erroneously as we now know, that the thermal conductivity remains 
finite in the limit of infinite system size. Simulations of the disordered latti ce con- 
nected to white noise reservoirs were carried out by Payton et al. |l53j |. They 
were restricted to small system sizes (N ~ 400) and also obtained a finite thermal 
conductivity. 

Possibly the first paper to notice anomalous transport was that by Matsuda 
and Ishii (MI) [13]. In an important work on the localization of normal modes 
in the disordered harmonic chain, MI showed that all high frequency modes were 
exponentially localized. However, for small u the localization length in an infinite 
sample was shown to vary as u>~ 2 , hence normal modes with frequency u> ~ u>d 
have localization length greater than N, and cannot be considered as localized. 
For a harmonic chain of length N, given the average mass m = {mi), the variance 
a 2 = {(mi — m,) 2 ) and interparticle spring constant k, it was shown that 



km \ 



1/2 



W < 80 > 



They also evaluated expressions for thermal conductivity of a finite disordered 
chain connected to two different bath models, namely: 

• model(a): white noise baths and 

• model(b): baths modeled by semi- infinite ordered harmonic chains (Rubin's 
model of bath). 

In the following we will also consider these two models of baths. For model(a) 
MI used fixed BC (boundary particles in external potential) and the limit of weak 
coupling to baths, while for case (b) they considered free BC (boundary particles 
not pinned) and this was treated using the Green-Kubo formalism given by Allen 
and Ford [l3|]. They found a = 1/2 in both conclusion which we will see 

is incorrect. The other two important theoretical papers on heat conduction in 
the disordered chain were those by Rubin and Greer [281 ] who considered model(b) 
and of Casher and Lebowitz [291 ] who used model(a) for baths. A lower bound 
[J] > 1/iV 1 / 2 was obtained for the disorder averaged current [J] in refn. [2cj ] who 
also gave numerical evidence for an exponent a = 1/2. This was later proved 



rigorously by Verheggen 3l(. On the other hand, for model(a), [29] found a rigorous 
bound [J] > 1/iV 3 / 2 and simulations by Rich and Visscher [651] with the same 
baths supported the corresponding exponent a = —1/2. The work in [32| gave a 
unified treatment of the problem of heat conduction in disordered harmonic chains 
connected to baths modeled by generalized Langevin equations and showed that 
models(a,b) were two special cases. An efficient numerical scheme was proposed 
and used to obtain the exponent a and it was established that a = —1/2 for 
model(a) (with fixed BC) and a = 1/2 for model(b) (with free BC). It was also 
pointed out that in general, a depended on the spectral properties of the baths. We 
will briefly describe this formulation [s^] here and see how one can understand 
the effect of boundary conditions on heat transport in the disordered chain. We 
will consider both the white noise [model(a)] and Rubin baths [model(b)]. One of 
the main conclusions will be that the difference in exponents obtained for these 
two cases arises from use of different boundary conditions, rather than because of 
differences in spectral properties of the baths. 
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The Hamiltonian of the mass-disordered chain is given by 

H = E it + S 2 MX/+1 " X ' )2 + 2*^ + ^ • (81) 
i=i ' i=i 

The random masses {tti-;} are chosen from, say, a uniform distribution between 
(m — A) to (m + A). The strength of onsite potentials at the boundaries is k! . 
The particles at two ends are connected to heat baths, at temperature Tl and Tr, 
and modelled by generalized Langevin equations. The steady state classical heat 
current through the chain is given by: 

J= k B (T L -T R) (82) 

where T N {u) = AT 2 (uj)\G 1n {uj)\ 2 , G(u) = Z~ 1 /k 
and Z = \-uj 2 M + £> - £(w)]/k , 



where M and <£ are respectively the mass and force matrices corresponding to 
Eq. (|8ip . As shown in [32|, the non-zero elements of the diagonal matrix for 
models(a,b) are En(w) = £attv(u;) = and given by 

= —i~fuj model(a) 
E(u) = k{l - mto 2 /2k - iu(m/k) 1/2 [l - muj 2 /(4/c)] 1/2 } model(b) , (83) 



where 7 is the coupling strength with the white noise baths, while in case of Rubin's 
baths it has been assumed that the Rubin bath has spring constant k and equal 
masses m. As noted above, T^{uj) is the transmission coefficient of phonons through 
the disordered chain. To extract the asymptotic N dependence of the disorder 
averaged current [J] one needs to determine the Green's function element Guv(w). 
It is convenient to write the matrix elements Z\\ = —m\oj 2 /k + 1 + k'/k — T,/k = 
—miio 2 /k+2—Tj' where £' = Yl/k—k' /k+1 and similarly Z^n = —mNU 2 /k+2—Y 1 l . 
Following the techniques used in [2^, [32|] one gets: 

\G 1N (to)\ 2 = Ar 2 |AjvMr 2 with (84) 
A n (lj) = D hN - E'CA.jv + £>i,jv-i) + S /2 D 2jiV -i 

where An(lv) is the determinant of Z and the matrix elements Di m are given by 
the following product of (2 x 2) random matrices ?}: 



where f, = ( 2 " ^f'" ^) 



We note that the information about bath properties and boundary conditions are 
now contained entirely in E'(o;) while D contains the system properties. It is known 
that l-Dj m | ~ e cNu) [27|], where c is a constant, and so we need to look only at the 

low frequency (uj ~ 1/iV 1 / 2 ) form of £' . Let us now discuss the various cases. 
For model(a) free BC correspond to k' = and so S' = 1 — ijuj/k while for 
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Figure 1. Plot of [J] versus N for free and fixed boundary conditions. Results are given for both 
models(a,b) of baths. The two straight lines correspond to the asymptotic expressions given in 
Eqs. J87I8SH and have slopes —1/2 and —3/2. Parameters used were m = 1, A = 0.5, fc = 1, 7 = 1, 

Tjj = 2, Tr = 1 and k a = 1 (from [33^ . 



model(b) free boundaries corresponds to k' = k and this gives, at low frequencies, 
£' = 1 — i{rnjk^l 2 u). Other choices of k' correspond to pinned boundary sites 
with an onsite potential k x 2 /2 where k Q = k' for model(a) and k Q = k' — k for 
model(b). The main difference, from the unpinned case, is that now i?e[S'] 7^ 1. 
The arguments of [32I ] then immediately give a = 1/2 for free BC and a = —1/2 for 
fixed BC for both bath models. In fact for the choice of parameters 7 = (m/c) 1 / 2 , 
the imaginary part of £' is the same for both baths and hence we expect, for large 
system sizes, the actual values of the current to be the same for both bath models. 
This can be seen in Fig. (pQ) where the system size dependence of the current for 
the various cases is shown. The current was evaluated numerically using Eq. (|82p 
and averaging over many realizations. Note that for free BC, the exponent a = 1/2 
settles to its asymptotic value at relatively small values (N ~ 10 3 ) while, with 
pinning, one needs to examine much longer chains (N ~ 10 5 ). 

These results clearly show that, for both models(a,b) of baths, the exponent a is 
the same and is controlled by the presence or absence of pinning at the boundaries. 
The reason that both models give the same exponents is that the imaginary part 
of their self energies, given by Eq. (183 j) . have the same small oj dependence. If the 
imaginary part of the self energy, i.e. T(o;), has a different u dependence, then 
one can get different exponents for the same boundary conditions [321 ] . 

For the present case, some more specific predictions can also be obtained. As 



mentione d be fore only modes u> ~ u>d are involved in conduction. An observation 



made in [1061 ] was that in this low frequency regime one can approximate [7jv(u;)] 



by the transmission coefficient of the ordered chain T®(uj). One can then write: 

[J]~(T L -T R ) T${u)du. 
Jo 

By looking at the iV — > 00 limit results for the ordered lattice, the following results 
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are obtained for the two different boundary conditions [331 ] : 

where c = 2^(mk) 1 / 2 / (j 2 + mk), 1 for model(a), model(b) respectively. For fixed 
boundaries we have c' = ^(mk) 1 ^ 2 /k Q 2 , mk/k 2 for model(a), model(b) respec- 
tively. A, A 1 are constant numbers, taken to be fitting parameters. For model(b) 
this agrees with an exact expression for [J}ft due to Papanicolau (apart from a fac- 
tor of 2tt) and with A = tt 3 / 2 / °° dt [t sinh(vrt)]/[(t 2 + 1/4) 1 / 2 cosh 2 (vrt)] « 1.08417 
(see [&]]). 

In the case where all sites of the chain are pinned (i.e. in the presence of a 



substrate potential) it was been noted in [32l . 1341 ] that the current decays exponen- 



tially with N and this was proved in [1561 ]. Another interesting result obtained in 
[33! ] is the case with a finite number n of pinned sites. It was shown, using heuristic 
arguments and numerics, that a = 3/2 — n for 2 < n « N. 

A question that has been discussed in the literature is that of uniqueness of 
the steady state. This depends on the choice of heat baths as well as the system 
studied. For models(a,b) of baths, the uniqueness of the steady state of a chain 
has been discussed in [28j and [291 ] . For baths consisting of harmonic oscillators 
one obvious necessary condition for uniqueness is that the bath spectrum should 
include the modes of the system (see for example [HI]). 



The quantum mechanical disordered chain has been discussed in [42j] where it 
was argued that the asymptotic system size dependence of the current should 
remain unchanged from the classical case (unlike the low temperature behaviour 
for ordered case). The temperature profile in a quantum mechanical disordered 
chain and quantum aspects such as entanglement have been numerically studied 
in 



3.4-2. Two dimensional disordered harmonic lattice 

So far there has not been much progress in understanding heat conduction in 
higher dimensional disordered lattices. The LEGF theory is still applicable and 
provides a general expression for the steady state current, Eq. (|64p . in terms of the 
phonon Green's function of the harmonic lattice. However in one dimension one 
could make progress by writing the transmission coefficient in terms of a product 
of random matrices as in Eqs. (j84|85j) . This enables one to use some known math- 
ematical theorems from which some analytic results could be obtained. Also this 
representation makes it possible to evaluate the current by very efficient numerical 
procedures. However in two and higher dimensions things become more compli- 
cated and it has not been possible to make much analytic progress. This state of 
things is also reflected in the fact that while in one dimension it can be proved 
exactly that all finite frequency states are localized, there is no such proof in two 
dimensions, although this is the general belief. As far as phonon localization is con- 
cerned a renormalization group study by John et al. [13] found that things are very 
similar to electronic localization. The important difference is again at low frequen- 
cies where one gets extended states. Their study indicates the following: in ID, all 
modes with to > 1/L 1 / 2 are localized; in 2D, all modes with uj > [log(L)] _1//2 are 
localized; in 3-D, there is a finite band of frequencies of non-localized states. How- 
ever this study was unable to extract the system size dependence of heat current 
in a disordered lattice in any dimension. 
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There have been a few simulation studies on the 2D disordered harmonic lat- 
tice. Lei Yang [591 ] considered a lattice with bond-missing defects and looked at 
the system size dependence of conductivity at various defect densities. The first 
observation that was made was that disorder gives rise to a temperature gradient 
across the system, unlike the flat profile for the ordered case. Also it was found 
that at small densities, the conductivity diverged logarithmically, while at larger 
densities, a finite heat conductivity was obtained. However this conclusion is prob- 
ably incorrect since the paper uses Nose-Hoover thermostats and it is know n tha t, 
for harmonic systems, these have equilibration problems 37, H, 60, 61, 6^, 155| |. 



The most detailed study of the 2D disordered harmonic lattice is by Lee and 
Dhar [g^]. They considered stochastic heat baths and looked at the case of mass 
disorder. In their model, the masses of exactly half the particles on randomly cho- 
sen sites of a L x L square lattice were set to one and the remaining to two. To see 
the effect of spectral properties of baths, two kinds of baths were studied: one with 
uncorrelated Gaussian noise and the other with exponentially correlated Gaussian 
noise. Simulations in disordered systems have to be done with care since one can 
have slow equilibration. In [g^] the authors first checked their simulation results 
by comparing them with results obtained from an exact numerical solution of the 
general RLL matrix equations Eq. (|45jl for a 8 x 8 lattice. The comparision, for 
the temperature at various lattice points, obtained by the two methods is shown 
in Fig. ([2]) and one can see excellent agreement between the exact results and 
simulation. Also shown in the inset are the disorder-averaged temperature profiles 
across the system for different sizes. One can see approach to a linear profile. To 
find the system size dependence of the current, lattices with sizes upto L = 256 
were studied. In Fig. ([3]) we show the data for the disorder averaged current for 
different system sizes and for the two different bath models. From this data the 
exponents a ~ 0.41 for white noise baths and a ~ 0.49 for the correlated bath 
was obtained. A special case of correlated disorder, first discussed in [30'], was also 
studied. Here the lattice was disordered in the conducting direction, but ordered in 
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Figure 3. Plot of disorder-averaged-current versus system size for a L X L lattice two different heat 
baths. The case of correlated disorder with white noise is also shown. For the full disorder cases, the solid 
lines are fits to the last three points and have slopes —0.59 and —0.51. For the case of correlated disorder, 

the slope from exact numerics (and also simulations) is compared to —1.5 which is what one expects 

analytically (from 1631 ). 



the transverse direction. Using the same methods as discussed in sec. ()3.3.2p . one 
can transform this into an effectively ID problem and then use the numerical tech- 
niques for evaluating the heat current as in the ID case. In [63-], the current was 
evaluated numerically (upto L = 512) as well as from simulations (upto L = 128 
and with white noise baths). Excellent agreement between the two confirmed the 
accuracy of the simulations. From the transformation to an effective ID problem 
it is possible to argue for an exponent a = —1/2 for the correlated disorder case 
and this could be already verified at the system size L = 512. This is somewhat 
surprising considering the fact that in the ID case, one has to go to sizes ~ 10 5 to 
see this exponent [see Fig. (JT])]. This result gives some confidence that the results 
obtained, for the fully disordered cases, are also close to the asymptotic values. An 
interesting observation in [g^] is that equilibration times for local temperatures is 
typically much larger than that for the current. This is expected since the tem- 
perature gets contributions from all modes including the localized ones which are 
weakly coupled to the reservoirs. On the other hand the current is mainly carried 
by low frequency extended modes. This point was also noted for the disordered ID 
chain in 1341. 



3.5. Harmonic lattices with self- consistent reservoirs 

As another application of the LEGF formalism we consider the problem of heat 
transport in a harmonic chain with each site connected to self-consistent heat 
reservoirs. The classical version of this model was first studied by Bolsteri, Rich and 



Visscher 64, 65] , who introduced the self-consistent reservoirs as a simple scattering 
mechanism for phonons which might ensure local equilibration and the validity of 
Fourier's law. The extra reservoirs connected to the system can roughly be thought 
of as other degrees of freedom with which the lattice interacts. It is interesting 
to note that the self-consistent reservoirs are very similar to the Buttiker probes 



68l . |69( ] which have been used to model inelastic scattering and phase decoherence 



November 19, 2008 



19:56 



Advances in Physics 



reva 



33 



in electron transport. In the electron case they lead to Ohm's law being satisfied 
just as in the harmonic chain the introduction of self-consistent reservoirs leads to 
Fourier's law being satisfied. In fact it has recently been shown how one can obtain 
both Ohm's law and Fourier's law in an electron model by using self-consistent 
reservoirs represented microscopically by noninteracting electron baths [701 ] . 

The ordered harmonic lattice with self-consistent reservoirs was solved exactly 
by Bonetto et al [67], in arbitrary dimensions, who proved local equilibration and 
validity of Fourier's law, and obtained an expression for the thermal conductivity 
of the system. They also showed that the temperature profile in the wire was linear. 
The quantum version of the problem was also studied by Visscher and Rich [6^] 
who analyzed the limiting case of weak coupling to the self-consistent reservoirs. 
We will here follow the LEGF approach as given in [44] to obtain results in the 
quantum-mechanical case. The classical results of Bonetto et al are obtained as the 
high temperature limit while the quantum mechanical results of Vischer and Rich 
are obtained in the weak coupling limit. We will consider only the one-dimensional 
case. A generalization to higher dimensions can be easily achieved, as in Sec. f|3.3|) . 

Consider the quantum-mechanical Hamiltonian: 



where we have chosen the boundary conditions xq = xn+i = 0. All the particles are 
connected to heat reservoirs which are taken to be Ohmic. The coupling strength 
to the reservoirs is controlled by the dissipation constant 7. The temperatures of 
the first and last reservoirs are fixed and taken to be T\ = Tl and T/v = Tr. For 
other particles, i.e , I = 2,3...(N — 1), the temperature of the attached reservoir 
Ti is fixed self-consistently in such a way that the net current flowing into any of 
these reservoirs vanishes. The Langevin equations of motion for the particles on 
the wire are: 

mxi = -mul(2xi - xi-i - xi + i) - muQXi - 7x1 + rji l = l, 2...N, (90) 

where the noise-noise correlation, from Eq. (|62p and with T(u>) = 70; for Ohmic 
baths, is given by: 

^( rn{v)Vm(u') + m(u')Vm(uj) ) = ^ coth( 2 ^J 5(uj + u') 5 im . (91) 

From the equations of motion it is clear that the I th particle is connected to a bath 
with a self energy matrix S^~(w) whose only non vanishing element is (S^)/; = ijuj. 
Generalizing Eq. (|64p to the case of multiple baths, one finds that the heat current 
from the Z th reservoir into the wire is given by: 




(89) 




(92) 



where T lm 
and G + 



4 Tr[ G + (uj) Ti(u)G'(u) F m (u)] 

[ -^M + d-^s+H]- 1 , 



= Jm[E+] . 



Here Ti m is the transmission coefficient of phonons from the I to the m reservoir. 
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Using the form of Ti one gets, in the linear response regime: 

Jl = 7 2 r cko ^ f; I [G+MW | 2 ( Tl - T m ) . (93) 

J — CO ^ i 

m=l 

For a long wire (N >> 1), for points far from the boundaries of the wire (I = yN 
where y = 0(1), 1 — y = 0(1)), one can explicitly evaluate the Green's function 
and show that: 



^ -a\l—m\ 

2muj 2 . sinh a ' ' ' 



where e a = z/2 ± [{z/2) 2 - l] 1 / 2 with z = 2 + lo^/lo 2 - oj 2 /uj 2 - ijlo / \muo c ) 2 , and 
we choose the root a such that i?e[a] > 0. Using this form of G~i m and assuming a 
linear temperature profile given by 

T l = T L + j^j(T R -T L ) , (95) 

one can see at once that, for any point I in the bulk of the wire, the zero-current 
condition J\ = is satisfied since X/m=-ooC — "2-)|e _a '' _m ' | 2 = 0. For points which 
are within distance 0(1) from the boundaries the temperature profile deviates from 
the linear form. Knowing the form of the temperature profile T; and the form of 
G^~ m , one can proceed to find the net current in the wire. It is easiest to evaluate the 
following quantity giving current Jij+i on the bond connecting sites I and (I + 1): 



' 1+1 = ^ {XlXl+l) = ~— U {^t) COS6Ch *2Jfe fl T 

N 

Y, k B T m Im{[G + {uj)\ lm [G + (uj)\* l+1 m } . 



N 

X 

m=l 



Using Eqs. (|94|95p one finally gets the following expression for the thermal con- 
ductivity k = JN/AT (obtained in the large N limit): 

-yk B f°° , to ( \ 2 , 2/ Hlj . / 1 1 \ . . 

k = — -At- / 9 cosech 2 (- ) — , (96) 

lomu>£m J-oo sinh or \2k B T J 2k B T \sinha sinha*/ 

where or(u;) = i?e[a], aj(uj) = Im[a]. In the high temperature limit 
{huj/2k B T) 2 cosech 2 (hiu/2k B T) -» 1, this gives, after a change of variables from 
w to q/, the following result for the classical thermal conductivity: 

2k B mto 2 (2 + v 2 ) rl 2 , sin 2 (aA 

Kd = / dai 



77r J (2 + v 2 ) 2 - 4 cos 2 (a/) 

k B muj 2 . 



7 (2 + v 2 + {v 2 (4 + v 2 )} 1 / 2 ) 



(97) 



where v = u>q/uj c . This agrees with the result obtained in refn. [63]. Another inter- 
esting limiting case is the case of weak coupling to the reservoirs (7 — ► 0). In this 
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case Eq. (|96j) gives: 

= fej Wo sm a/ cosech ( 2^r } ' (98) 

where w a = Wg + 2w 2 [l — cos(a/)] . 

This agrees with the result obtained in [6a ]. In the low temperature limit, Eq. (|98D 
gives k« ~ e -t^ /k B T jrpi/2 for Wo _^ q and ^ T for Wo = q_ Ag notec } in ]g 

the expression for thermal conductivity (in the weak scattering limit) is consistent 
with a simple relaxation-time form for the thermal conductivity. The temperature 
dependence of n wc then simply follows the temperature dependence of the specific 
heat of the ID chain. 

In the general case where the coupling constant has a finite value, the low- 
temperature behaviour again depends on whether or not there is an onsite poten- 
tial. The form of the low-temperature behaviour is very different from the case of 
weak coupling. For small T it is easy to pull out the temperature dependence of 
the integral in Eq. (|96p and one finds that k ~ T 3 for v 7^ and k ~ T 1 / 2 for 
v = 0. 

A nice extension of the above problem has been done by Roy pFlj ] who considered 
the case where the coupling of the bulk lattice points to reservoirs (7') was taken to 
be different from that of the boundary points to the end reservoirs (7). The quan- 
tum mechanical case with Ohmic reservoirs and the linear response regime were 
studied. For small values of 7' the transition from ballistic to diffusive transport 
could be seen with increasing system size. It was shown that the current, for any 
system size, could be written in the form 

k(T)AT , . 

where £ was an effective mean free path for phonons which depended on 7'. Thus for 
N « £ one gets ballistic transport while for N >> £ one gets diffusive transport. 
In Fig. (]4|) the temperature profiles for different system sizes is shown and one can 
see the transition, from a relatively flat (in the ballistic regime), to a linear profile 
(in the diffusive regime). 



4. Interacting systems in one dimension 

In the case of systems with interactions there are few analytic results, for one- 
dimensional interacting particle systems, and these are all based on use of the 
Green-Kubo formula. All these theories aim at calculating the equilibrium current- 
current correlation function C(t) = (J(0)J(t)). Mainly there are three different 
theoretical approaches: renormalization group theory of hydrodynamic equations, 
mode coupling theory and the Peierls-Boltzmann kinetic theory approach. We will 
discuss these in Sec. (I4.ip . We will also discuss some exact results, which have been 
obtained for certain models for which the dynamics is not completely deterministic 
but includes some stochastic component. All the analytic approaches predict that 
momentum conserving systems in ID exhibit anomalous transport with conductiv- 
ity diverging as a power law k ~ N a . However there is disagreement on the precise 
value of a and the number of universality classes. We will present the results of 
simulations for various momentum conserving models in Sec. (|4.2.1j) . In general 
one finds anomalous transport with k ~ N a and again there is disagreement in 
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Figure 4. Plot of the temperature profile {T; } as a function of scaled length l/N for different N with 
7 = 1.0 and 7' = 0.1. Here mean free path I ~ 30 (from |71j1. 



the value of a obtained in simulations by different groups. However, the evidence 
seems strong that there is a single universality class with a = 1/3. 

For momentum non-conserving interacting systems the prediction of theory is 
that Fourier's law is valid. Simulations also show [Sec. (|4.2.2p ] that the presence 
of interactions (nonintegrable) and an external substrate potential are sufficient 
conditions to give rise to a finite thermal conductivity. Note that momentum non- 
conserving models are a bit unphysical in the context of heat transport by molecular 
motion except when one is considering wires and thin films attached to substrates. 
In the case of electron transport the background ionic lattice naturally provides 
an external potential for the electron and so momentum non-conservation makes 
physical sense. We will also briefly discuss a class of rotor models that have been 
studied in simul ations wh ich show Fourier behaviour inspite of absence of an ex- 
ternal potential 1351 . 13d ] . These models are somewhat special, in the sense that it 
is more natural to think of the phase space degrees of freedom as local angle vari- 
ables and thus these models should probably be thought of as angular-momentum 
conserving rather than linear-momentum conserving models. 



4.1. Analytic results 

4- 1.1. Hydrodynamic equations and renormalization group theory 

This approach was first proposed by Narayan and Ramaswamy [72]]. Here one 
first notes that a one dimensional (ID) system of interacting particles will, at 
sufficiently large length scales, behave like a fluid. Suppose that the only conserved 
quantities in the system are the total number of particles, the total momentum and 
total energy. One can then write hydrodynamic equations to describe the variation, 
in time and space, of the density fields corresponding to these conserved quantities. 
Namely we have p(x, t), g(x, t) = p(x, t)v(x, t), and e(x, t) for number, momentum 
and energy densities respectively and where v(x, t) is the local average velocity 
field. This basically gives the Navier-Stokes equations for a ID fluid. After adding 
noise terms to account for thermal fluctuations in the system the equations are 
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given by 

d t p + d x (pv) = 0, 
d t (pv) + d x (pv 2 ) = -d x p + (d 2 v + rj v , 

dte + d x [(e + p)v] = k8 2 x T + ([(d x v) 2 + vd 2 x v] + Ve , (100) 

where the noise terms satisfy {r] v (xi,ti)r] v (x' 1 ,t' 1 )) oc —ksT5(xi — Xi)S(ti — t[) 
and similarly for r] 6 . The local temperature T(x,t) and pressure p(x,t) are im- 
plicit functions of p and e. There are two transport coefficients, the viscosity £ 
and the thermal conductivity k. The above equations can be solved in the linear 
approximation and this gives C(t) ~ t~ l l 2 which in turn implies a divergence of 
the conductivity. However this divergence also means that the linear approxima- 
tion is not good and one needs to take into account the nonlinear terms in the 
Navier-Stokes equations in order to get the correct long time behaviour of various 
correlation functions. In general this would require a RG analysis but it is argued 
in [73] that the exponents can be obtained from symmetry considerations. Using 
the Galilean invariance of the system and the fact that equal time correlations obey 
equilibrium statistical mechanics they finally obtain: 

C(t) ~ t~ 2 / 3 . (101) 

For a finite size system, using the arguments in Sec. (|2.ip . one puts a upper cut-off 
~ N in the Green-Kubo integral and this then gives k ~ iV 1 / 3 . Thus a = 1/3. 
Note that in this treatment, the details of the form of the Hamiltonian are unim- 
portant. The only requirements are the presence of the three conservation laws and 
also, the interactions should be such that the nonequilibrium state satisfies local 
thermal equilibrium and should be describable by coarse grained hydrodynamic 
equations. We note that the possibility of breakdown of hydrodynamic equations 
in a one-dimensional fluid system has recently been pointed out [t3 ] . 

An interesting question that arises in the context of the hydrodynamic theory 
is the behaviour of the other transport coefficient in the equations, i.e. , the bulk 
viscosity £. This has not been investigated much except in the work in refn. [75|] 
who, somewhat surprisingly, find that this transport coefficient is finite. 

4- 1.2. Mode coupling theory 

This approach was first applied in the context of heat conduction by Lepri, 
Livi and Politi 18211 and has subsequently been used by several other authors [82j, 
[H, l84l. l85l . [H, H3, [H, H^]. We will here outline the main steps as described in 
refn. [87]]. Mode coupling theory (MCT) again begins with the realization that 
the divergence of conductivity is a result of the long time tails of the current- 
current correlation function which in turn can be attributed to the slow relaxation 
of spontaneous fluctuations of long-wavelength modes in low dimensional systems. 
For a ID oscillator chain with periodic boundary conditions one considers the 
normal mode coordinates of the harmonic lattice which are given by: 

1 N 

Q{q) = -7= ^2 Xn ex P(-^) , (102) 

* n=l 



where the wavenumber q = 27rk/N with k = —N/2 + 1, —N/2 + 2, ...N/2 (for even 
N). The evolution of a fluctuation at wavenumber q excited at t = is described 
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by the following correlation function: 



_ (Q*( g ,t)Q(g,0)) 
{q,) ~ (\Q(Q)\ 2 ) ' ( } 

The long time decay rate of this quantity at small q is the main object of interest 
and MCT is one approach to obtain this. Basically one writes a set of approximate 
equations for G(q,t) and this is then solved self-consistently. Formally one can in 
fact write an exact equation for the time evolution of G(q,t). Using the Mori- 
Zwanzig projection methods one gets the following equation: 



G(q, t)+£ f F(q, t - s)G{q, s) ds + LU 2 (q)G(q, t) = , (104) 
J 

where the memory kernel T(q, t) is proportional to (J~(q, t)J-(q, 0)), with J-(q) being 
the nonlinear part of the fluctuating force between particles. The coupling constant 
e and the frequency u>(q) are temperature dependent input parameters which have 
to be computed independently. Equations (|104p must be solved with the initial 
conditions G(q, 0) = 1 and G(q, 0) = 0. The mode-coupling approach proceeds by 
replacing the exact memory function T(q, t) with an approximate one, where higher 
order correlators are written in terms of G(q, t). Consider now the FPU interaction 
potential U(x) = k^x 1 /2 + ksx^/3 + k^x^/A. In the generic case, in which k^ is 
different from zero, the lowest-order mode coupling approximation of the memory 
kernel gives: 

T(q,t)= u; 2 (q)^ Yl G(p,t)G(p',t) . (105) 

P+P'— <7=0,±7T 



On the other hand for the case k% = 0, k^ ^ 0, one gets: 



T(q,t) = u 2 (q)(^) 2 E G(p,t)G(p',t)G(p",t) . (106) 



-2tt\2 

P+p'+p"— g=0, ±7T 



Here p ,p',p" range over the whole Brillouin zone (— vr, tt). Using either of Eq. (|105p 
or Eq. (|106j) in Eq. (I104p gives a closed system of nonlinear integro-differential 
equations. The coupling constant e and the frequency cu(q) are taken as parameters 
which can be obtained from the harmonic approximation. The solution of these 
equations again involves making a number of other approximations and the final 
result one obtains for G(q, t) at small values of q is the following form: 

G(q,t) = A(q,t)e^ q)t + c.c 

where « q .t) - l*?**? b ' k ^° (107) 
\g{e l/2 tq 2 ) for A; 3 = 0, k 4 + 0. V ' 

Finally one can relate the current-current correlation function C(t) to the correlator 
G(q,t). Again making the same approximation of retaining only the lowest order 
correlation functions one gets: 



C(t)«^G 2 (9,i) , 
q 



(108) 
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and plugging into this the result from Eq. (I107p . one finally obtains: 

rm J*"' 73 for hm\ 

C V~\t-V> for k 3 = 0,^0. (1 ° 9) 

Inserting this into the Green-Kubo formula (with a cut-off proportional to N) then 
gives us k ~ A^ 1 / 3 for the odd potential and k ~ N 1 ! 2 for the even potential. 

Another recent MCT study is by Wang and Li [H, d^] who look at the effect of 
transverse degrees of freedom on the value of a. They consider a one dimensional 
chain where the particle positions are now two-dimensional vectors instead of being 
scalar variables. The Hamiltonian they consider corresponds to a polymer with 
bending rigidity and is given by: 

H = E || + it (i r '+i - r <i - a ) 2 + ^ c ° s (<w > ( n °) 

where {17, p^}, for I = 1,2...^ denote two dimensional vectors and cos(0/) = 
— n;_i.ri2 with n; = Ar;/|Ar;| and Ar; = ri + \ — 17. Based on their MCT anal- 
ysis they suggest that the generic effect of including transverse degrees is to give 
a = 1/3, while for the purely longitudinal model one has a = 2/5. 

For momentum non-conserving systems MCT predicts a finite conductivity. 

4-1.3. Kinetic and Peierls-Boltzmann theory 

In the kinetic theory picture, one thinks of a gas of weakly interacting particles, 
which are the heat carriers. These heat carriers could be molcules in a gas, elec- 
trons in a metal or phonons in a crystal. Using the idea that the heat carriers are 
experiencing random collisions, and hence moving diffusively, one can do a simple 
minded calculation. This gives us a simple expression for the thermal conductivity, 
namely k ~ cvl, where c is the specific heat capacity per unit volume, v the typical 
particle velocity and t the mean free path of the particles between collisions. 

The Boltzmann equation approach gives a more systematic derivation of the re- 
sults of kinetic theory, and was first developed for the case of molecular gases. In 
this theory, one writes an equation of motion for the distribution function /(x, p, t), 
where /(x, p, t)d?xd?p (in 3D) gives the number of particles in the volume d 3 xd s p. 
The presence of collisions makes the Boltzmann equation equation nonlinear, and 
then one has to solve the equation under various approximations. The final result is 
quite often in the form of the kinetic theory answer, with an explicit expression for 
the mean free path £. For phonons, the Boltzmann theory of conductivity was de- 
veloped by Peierls [91(. He wrote the Boltzmann transport equation for the phonon 
gas and pointed out the importance of lattice momentum non-conserving processes 
(Umpklapp processes) in giving rise to finite conductivity. Solving the Boltzmann 
equation in the relaxation time approximation gives a simple kinetic theory like 
expression for the thermal conductivity, k ~ f (iqCqfqTq, where r q is the time 
between collisions, and q refers to different phonon modes of the crystal. The relax- 
ation time r q can get contributions from various sources, such as phonon-phonon 
interactions and impurity scattering, and its calculation from first principles is one 
of the main tasks. In three dimensional solids, the Peierls-Boltzmann theory is 



well-developed 92j, l93f] and probably quite accurate. One worry here is that the 
meaning of the distribution /(x, q, t) for phonons is not really clear, since phonons 
are extended objects. The recent work of Spohn [94j] tries to give a rigorous basis 
for the phonon Boltzmann equation for a crystal with a weakly anharmonic onsite 
potential. We note that, as far as making definite predictions (starting from a given 
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Hamiltonian) on the actual conductivity of a system and properties such as the 
temperature dependence of the conductivity, the kinetic theory approach probably 
has more chance of success than the other two approaches described before. 

Let us now consider the application of the kinetic theory approach to the one 
dimensional case. If we look at the result k ~ f dqc q v 2 T q , we see that a divergence 
with system size can arise because the relaxation time r q , and correspondingly the 
mean free path l q = v q T q , becomes large at small q. Let us assume that c q ,v q are 
constants, and that T q has the power-law dependence r q ~ q~ a . Then all modes 
with q < q m ~ L -1 / a travel ballistically and this immediately gives K ~ L 1_1 / a , 
and for a > 1 one would get a diverging conductivity. 

In a way the kinetic approach is similar to the MCT method. Here too, one tries 
to calculate the rate of decay of long-wavelength fluctuations at small q, but now 
using a different approximation scheme. The current-current correlation function is 
then related to this decay constant. The Kinetic theory approach for the FPU chain 
was first considered by Pereverzev [95[ who studied the model with k% = and a 
small non-zero value for k^. It was noted that the approximate time evolution of a 
fluctuation in the average energy e q of a mode with wavenumber q is given by the 
homogeneous classical linearized Peierls equation. This equation is then brought 
to the following form, corresponding to the relaxation time approximation: 

d{e g (t)) 1 

— "77 — = ((e?) - k B T) , (111) 

OX T q 

where one has an explicit form for r q . For small q, making some more approxima- 
tions enables one to evaluate r q and one finds r q ~ g~ 5 / 3 . Finally using the same 
set of approximations and in the limit N — > oo one can show that: 

nh.2 rp2 fir 

C (t) = \ dqe^vl , (112) 

7T Jo 

where v q is the phonon group velocity. At small q the phonon velocity v q ~ const 
and the above equation gives: 

C(t) ~ i~ 3 / 5 . (113) 

This then implies k ~ N 2 ^ . In fact, the arguments given at the beginning of this 
paragraph directly give this (putting a = 5/3), and one does not need to find C{t). 

The kinetic theory approach has been made more rigorous by the work of Lukkari- 
nen and Spohn [96j] . They also work with the linearized collision operator and make 
the relaxation time approximation, and for the quartic FPU chain they confirm the 
result in [95], namely C(t) ~ t~ 3 / 5 . However they point out the possibility that 
the kinetic theory approach may not be able to predict the correct long-time decay 
of the correlation function. Another paper using the linearized Peierls-Boltzmann 
equation for the quartic Hamiltonian also finds k ~ N 2 / 5 [97]. Finally a quantum 



calculation of the phonon relaxation rate at small q has been carried out in 98l.l9£ 
They studied both the cubic and quartic FPU chain and obtained relaxation times 
" 3 / 2 , g -5 / 3 for the two cases respectively. 



For the case of mo mentu m non-conserving systems, Lefevere and Schenkel 13Cll | 



and later Aoki et al. [l31 ] have used the kinetic theory approach for the case of 



weak anharmonicity and obtained a finite conductivity. 
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4- 1-4- Exactly solvable model 

A harmonic chain with a energy conserving stochastic dynamics was considered 
by Kipnis et al. [7a] who could prove exactly that the model satisfied Fourier's law. 
The dynamics was momentum non-conserving and completely stochastic so it is not 
surprising that Fourier's law was obtained. Models with self-consistent reservoirs 
can also be viewed as stochastic models (but with Hamiltonian components in the 
dynamics) where energy is conserved on average, while momentum is not and again 
Fourier's law is satisfied. 

Recently a similar stochastic model, but in which total momentum conservation 



was also enforced, was introduced by Basile et al. [77], |7g|. In their lattice model 
the dynamics consisted of two parts. Apart from a deterministic Hamiltonian dy- 
namics the system was subjected to a stochastic dynamics which conserved both 
total energy and momentum exactly. The stochastic dynamics consisted of a ran- 
dom exchange of momentum between three neighboring particles (in ID) while 
conserving both energy and momentum. Thus a triplet of particles with momenta 
(Pi— liPliPl+\) is chosen and this set performs a diffusive motion on the curve given 
by: 

Pl-1 +Pl +P1+1 = const. 

Pl-1 j_ P 2 l j_ Pl+l 

1 — H — = const. 

2m 2m 2m 

The Hamiltonian of the system was taken to be that of a harmonic system. A 
Fokker-Planck equation for the probability density P(x, p, t) could be written, 
which in ID is given by: 

BP 

— = (L H + 1 S)P, (114) 

where L H is the usual Liouville operator for the given Hamiltonian and S, the 
generator of the stochastic perturbation of strength 7, has the form 



S 



6 

I 




with Yi = (pi - Pi+i)d Pl _ t + (pi +1 -pi-i)d Pl + (pj_i -pi)d Pl+1 . 

The authors were able to compute exactly an explicit form for the current-current 
correlation function C(t), for system size N — > 00, and from this they found the 
following asymptotic long-time behaviour: 

. . for no pinning 

C W~^-3/2 - + i • ■ ( 115 ) 

with pinning. 

Plugging this into the Green-Kubo formula one gets a = 1/2 in the unpinned case, 
while for the pinned case a finite conductivity is obtained. One can argue that the 
stochastic dynamics in a way mimics anharmonicity and the problem considered 
corresponds to an even interaction potential. The latest prediction from MCT also 
gives a = 1/2, which agrees with the result from this model. However all simulation 
results of momentum conserving interacting Hamiltonian models give exponents 
quite far from this value (between 0.3 — 0.4). The exponent a = 1/2 then comes as 
quite a surprise. One possibility is that the choice of a harmonic Hamiltonian makes 
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the model special, and thus solvable, and at the same time makes it a non-generic 
case. 

Some simulations with a dynamics which is roughly similar to the above stochas- 
tic dynamics were recently done with FPU type anharmonic terms included in the 
Hamiltonian [7§]. These are equilibrium simulations using the Green-Kubo for- 
mula. The authors have argued that their results support the two-universality 
class scenario. It will be interesting to understand in more details the role of the 
Hamiltonian part of the dynamics in determining the exponent a in this model. 

More recently, a similar one-dimensional stochastic model with random two par- 



ticle momentum exchanges has been numerically [80] and analytically [8l|] studied, 
for the nonequilibrium case with Langevin heat baths attached at the two ends. 
Apart from confirming the exponent a = 1/2, these studies have also looked at 
the temperature profiles. An analytic expression for the temperature profile was 
obtained and it was noted that the profiles were very similar to those obtained for 
FPU chains. 



4.2. Results from simulation 

4-2.1. Momentum conserving models 

Gas of elastically colliding particles of two masses: One of the simplest 
model of interacting particles that one can consider is a gas of elastically collid- 
ing point particles where the boundary particles interact with thermal reservoirs, 
usually modeled by Maxwell boundary conditions. If all the particles h aye e qual 
masses then this model, without reservoirs, is the so-called Jepsen model As 



far as heat conduction properties are concerned the model is somewhat trivial. This 
is because at each collision the particles simply exchange momentum and so the 
net heat transfer can be calculated by considering a single particle that is bouncing 
between the hot and a cold walls. One finds a system- size independent heat cur- 
rent J = 4 /2 (2m/7r) 1 /2p(T2r R -r|T L )/(rJ /2 T R + r^ /2 r L ), where p is the particle 
density, and a flat temperature profile given by T = (T^Tr) 1 / 2 . Thus this model is 
somewhat like the ordered harmonic chain. However the model becomes interesting 
and non-trivial if one considers a dimerized model where alternate particles have 
different masses say m\ and mi. In this case one finds a current which decays with 
system size, and a slowly varying temperature profile. 



The diatomic hard particle gas model was first studied by Casati 10l| but the 
numerical results were not sufficient to draw any definite conclusions. T his m odel, 
along with the diatomic Toda lattice, were later studied by Hatano |l02j |. Us- 
ing nonequilibrium simulations and system sizes upto N = 5000, an exponent 
a rj 0.35 was obtained for both these models. The current-current correlation 
function was also evaluated for a periodic closed system and it was found that 
C(t) ~ j\f-o.65 cons i s tent with the nonequilibrium results. Subsequently, a num- 
ber of further studies were made using both nonequilibrium simulations, and also 
the Kubo formalism and using much larger system sizes. Unfortunately there is 
not much agreement on the numerically obtained value of the exponent. The 
vario us reported va lues include: Garrido et al. (a = implying Fourier behav- 



ior) |Iq3, [Toi, [l05|], Dhar (a ra 0.2) [loi], Grassberger et al. (a *j 0.33) [l08l | 



bavin et al. (a ra 0.2) [Io3] and Casati et al. (a ~ 0.25) |ioi |. However, based on 



the theoretical predictions, there seems reason to believe that the value obtained 
by Grassberger et al. is the correct one and here we will di scuss their results in 



some detail. We also mention here the work of Cipriani et al. [llOl ] who performed 



zero-temperature studies on diffusion of localized pulses and using a Levy walk 
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Figure 5. (a): Log-log plot of J/(T<2 — Ti) versus N for four values of the mass ratio A. (b): Part of the 
same data divided by N a with a = 0.32, so the y-axis is much expanded (from 108]). 



W^. Apart 



interpretation concluded that a = 0.333 ± 0.004. 

We now present some of the results obtained by Grassberger et al. 
from looking at much larger system sizes (upto ./V = 16384), they made the obser- 
vation that the asymptotic behaviour is easier to observe at some optimal value of 
the mass ratio A = m^/fni. It was argued that A = 1 and A = oo were special 
integrable limits where one would clearly get ballistic and non-typical behavior. 
If the value of A was too close to 1 or too large then one would have to go to 
very large system sizes to see the correct asymptotic form. However by choosing 
an appropriate value of A, one can reach asymptotic behaviour much faster. This 
feature can be seen in Fig. ([5]) where the system size dependence of the current for 
different values of A is given. One can see that for A = 2.62, asymptotics is reached 
faster than for A = 1.618 and A = 5.0. The value of the exponent obtained from 
this data was a = 0.32^'^. Equilibrium simulations were also performed and in 
Fig. (JH) results are shown for the current-current autocorrelation function obtained 
for a periodic system. For large system sizes one can see a t -0 - 66 decay with a cutoff 
at t oc N . This again gives a = 0.34 in agreement with the value obtained from 
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Figure 6. Total heat current autocorrelation, i ' 66 ^- 1 ( J(t)J (0)) lor A = 2.2 and T = 2. Total 

momentum is P = (from [10811 ). 




the nonequilibrium simulation. 

Some interesting features were seen in the temperature profiles also and we now 
discuss these. Simple kinetic theory predicts that the thermal conductivity of a 
hard particle gas should have a temperature dependence k ~ T 1 / 2 (this can also be 
obtained from the Green-Kubo formula). Now if we plug this in Fourier's law then, 
with specified boundary temperatures, one easily obtains the following nonlinear 
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Figure 8. Log-log plot of the conductivity as a function of the number of particles for the random 
collision model introduced in [lllj |. The upper plot is for all the particles with mass ratio A = 1, whil e 
the lower plot is for A = 2.62. The slopes for the two plots are 0.29 ± 0.01 and 0.335 ± 0.01 (from HTllp . 



temperature profile: 

T k (x) = [T 3 L /2 (1 - x/N) + T^x/N] 2 / 3 . (116) 



In refn. 106fl . it was noted that a convergence, of the actual nonequilibrium tem- 



perature profile T(x), to the kinet ic th eory prediction T k (x) given above, seemed to 
take place. However the study in 1081 ] found that this apparent convergence stops 



as one looks at larger systems and in fact one finds T{x) — T k {x) attains a non-zero 
profile for N — > oo. This is shown in Fig. ([7]). This result indicates that there is a 
problem in writing Fourier's law in the form J = — katVT, with kn defined as a 
length dependent conductivity. We will see similar problems with other ID models. 

Random collision model: In studies on heat conduction, one often finds that 
a faster convergence to the asymptotic system size limit can be obtained by partic- 
ular choices of model and parameter values. In this c onte xt one should mention a 



stochastic model introduced by Deutsch and Narayan lll| . They consider a binary 
mass gas of hard point particles where any particle's motion is strictly confined to 
one dimensions while its momentum is a two dimensional vector (p x ,Py)- During 
a collision between two particles their momenta gets changed randomly while con- 
serving total energy and both components of total momentum. Physically one can 
think of this model as approximating a system of small particles with rough sur- 
faces, moving in a narrow tube. The results of nonequilibrium simulations for two 
different mass ratios is shown in Fig. ((H). For the case mxjmi = 2.62, at system 
sizes as small as N ~ 10 3 , one already gets a = 0.335 ± 0.01, which is close to the 
expected value a = 1/3. 

Fermi-Pasta-Ulam chain: The Fermi-Pasta-Ulam (FPU) model consists of an 
oscillator chain with harmonic as well as anharmonic nearest neighbour interparti- 
cle i nteractions. This model was first studied by the authors in a landmark paper 
112f | where they wanted to verify the common assumption of statistical mechan- 



ics that anharmonic interactions should lead to equilibration and equipartition. 
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Surprisingly the numerical experiments on the FPU chain gave a negative result, 
i.e. the chain failed to equilibrate. Understandin g th e FPU results led to the de- 
velopment of new areas and concepts in physics 113| j. It is probably fair to state 
that a complete understanding of the problem is still lacking. For example it is be- 
lieved that for high enough energy densities and large system sizes one may achi eve 
equilibration but there are details which are not yet understood precisely |ll4j |. 

What about heat transport across the FPU chain ? Clearly this is the simplest 
model to study in order to see the effect of anharmonicity on heat transport. 
One might suspect that the equilibration problem of the FPU chain is likely to 
show up in some way when one looks at transport properties, especially so when 
one thinks in terms of the Green-Kubo formula. It turns out that a FPU chain 
connected to heat reservoirs is better-behaved. It can in fact be proved rigorously 
that a FPU chain connected to equal temperature Langevin heat baths at its 
two ends will always equilibrate. At long times it will converge uniquely to the 
appropriate Boltzmann-Gibbs distribution. Also it can be shown that even with 
unequal temperature baths (s toch astic or Hamiltonian) the system reaches a unique 
nonequilibrium steady state jl!5l ] and this is very reassuring when one begins to 
consider heat transport studies in the FPU chain. 

The first study of heat conduction in the FPU chain was by Lepri et al. [lid ] 
who considered an interparticle potential of the form U(x) = k^x 2 jl + k^x^/A and 
performed nonequilibrium simulations with Nose-Hoover baths. Looking at system 
sizes upto N = 400 they obtained a = 0.55 ± 0.05. In a subsequent paper 
by studying systems upto N = 2048, they obtained a « 0.37. They also found a 
highly nonlinear and singular temperature profile and noted that this was true even 
for relatively small temperature differences applied to the ends. We will comment 
more on the temperature profil e of t he FPU chain later in this section. 

Since the important work of |lld |. a large amount of numerical and analytical 
work has been carried out on heat conduction in the FPU chain. We first summarize 
the various analytic results discussed in Sec. (14. ip . We assume that the interparticle 
interaction is of the general FPU form U(x) = k2X 2 /2 + /c3X 3 /3 + k^x 41 /^. The 
predictions from theory are then: 

(i) Renormalization group theory of hydrodynamic equations: This predicts that 
there is only one universality class with a = 1/3. 

(ii) Mode-coupling theory: This predicts that there are two universality classes 
depending on the parity of the leading nonlinearity in the Hamiltonian. For the 
case where the leading nonlinearity is cubic, i.e. k% ^ 0, the prediction is a = 1/3 
while for k% = 0, ^ / 0, the prediction is a = 1/2. 

(iii) Kinetic theory and the Peierls-Boltzmann equation approach: This gives 
q = 2/5 for the quartic case. 

Results o f simulat i ons: As we have seen in the last section simulations of hard par- 
ticle gases |l08l . 110 , 111 ] seem to indicate a value a = 1/3 for t he heat co nduction 



exponent, though even here the issue is not completely settled 107 . 109]. On the 



other hand, numerical simulati ons o f oscillator chains, including FPU chains, give 
various exponents @, [82, 84, 88, 1161 ] for different systems, often slightly higher than 



1/3. This seems consistent with early results from mode-coupling theory (MCT) 



which predicted a = 2/5 . The most recent MCT analysis [861 . 87 ] predicts that 
a = 1/2 for potentials U{x) with quartic leading nonlinearity while for potentials 
with cubic nonlinearity, there seems to be agreement between different theories 
about a = 1/3. Here we will focus on simulations for the even potenti al on ly. We 
will discuss the results of the most recent simulations by Mai et al. ll 18tl of the 
even potential FPU model and another simulation by Dhar and Saito [1571 ] of the 
alternate mass FPU chain. 
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Figure 9. Kinetic temperature profile for a FPU-/3 chain with N = 16384, T L = 2.0, T R = 0.5 and 
7 = 2.0. The temperature evaluated from the first three even moments of the velocity are shown; their 
agreement indicates a Gaussian velocity distribution. ( Inse t) Normalized temperature profiles for 

different N (from [HI). 



An important aspect addressed in the simulations by Mai et al. was the effect of 
boundary conditions. It is well known 0] that the of coupling of a system to thermal 
reservoirs leads to so-called contact resistances. These show up, for example, in the 
jumps that one observes in the temperature profile in such a system. It is only for 
sufficiently large system sizes, when the resistance of the system is much larger than 
the contact resistance, that one can neglect the contact resistance. In simulations 
where one is interested in determining the precise dependence of current on system 
size, it is important to ensure that one has reached the required system size where 
contact resistances are negligible compared to the actual syst em r esistance. This 
point has been discussed in some detail by Aoki and Kusnezov [117| . The study by 
Mai et al. ensured this by performing simulations with two different baths, namely, 
stochastic white noise baths and the deterministic Nose- Hoover bath. Further they 
did simulations for different coupling strengths of the system to reservoir. It was 
found that for small system sizes the current values were significantly different 
(for the same applied temperature difference). This is expected since the contact 
resistance, which is different for the different boundary conditions, dominates the 
transport current. However at large system sizes, the actual values of the currents 
for all the different cases tend to converge. In this system-size regime one is thus 
assured that boundary effects have been eliminated and one can then extract the 
correct exponent. 

The simulation in |l!8l | was done for parameter values hi = 1, = 0, = 1 and 
m = 1. The temperature at the two ends were fixed at Tl = 2.0 and Tr = 0.5. 
Both white noise baths with coupling parameter 7 and Nose-Hoover baths, with 
coupling parameter were studie d. Th e white noise simulations were done using 
a velocity- Verlet type algorithm 119l |. while the Nose-Hoover simulations were 



implemented using a fourth order Runga-Kutta integrator. Time steps of order 
dt = 0.0025 — 0.005 were used and, for the largest system size (N = 65536), upto 
10 9 equilibration steps and an equal number of data-collecting steps were used. 

The temperature profile for a chain of size N = 16384 is shown in Fig. ([9]). 
The temperature is defined through the first three even moments of the velocity 
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Figure 10. (Top) Conductivity versus N. The last five points fit to a slope of a = 0.333 ± 0.004. 
(Bottom) JN 1—a versus N for a = 1/3 and a = 2/5. In the large N regime, a is definitely less than 2/5 
and appears to agree quite well with the 1/3 prediction for all data sets. Langevin bath s with 7 = 0.4, 2 
and 10, and one data set with Nose-Hoover baths, are shown (from [11810 . 



as Ti = (tf), T x = ((uf)/3) 1 / 2 and T\ = ((vf)/2 - (vf) /SO) 1 / 3 respectively. Their 
agreement indicates that local thermal equilibrium has been achieved and the local 
velocity distribution is close to Gaussian. Also we notice that the boundary jumps 
are almost absent for this system size. The inset shows smaller system sizes where 
the boundary jumps, arising from the contact resistance, can be clearly seen. As 
noted and discussed earlier by Lepri et al.0| the temperature profile is nonlinear 
and this feature seems to be present even for small temperature differences and is 
another indication of anomalous transport. As for the hard particle case this also 
indicates that one cannot find the temperature profile using a temperature (and 
system size) dependent conductivity in Fourier's law. 

In Fig. (fTUl) (upper figure), the conductivity defined as n(N) = JN/ AT is plotted 
against system size. This data gives 

a = 0.333 ± 0.004 . (117) 

The results of various simulation runs with Langevin baths with different damping 
constants 7 = 0.4, 2, and 10 as well as the deterministic Nose-Hoover thermostat 
is shown (lower figure) in Fig JlOh . This compares the RG prediction (a = 1/3) 
and the old MCT prediction (a = 2/5) for systems with these different baths 
and bath parameters. As can be seen in the figure, an asymptotic exponent of 
1/3 is attained for all these systems, whereas the apparent exponents for smaller 
N depend on system parameters. It is possible to understand the deviation of 
the apparent exponent from 1/3 for small system sizes. As shown in refn. [2|, if 
the damping constant for the Langevin baths is very large or small, there is a 
large 'contact resistance' at the boundaries of the chain. The current only depends 
weakly on N, resulting in an apparent a > 1/3 (similar considerations apply to 
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Figure 12. Plot of the heat current J versus system size in the alternate m ass F PU chain for different 
values of the mass ratio A = 1.0, 1.1, 1.22 and 1.5 (from [157l | 1. 



Nose-Hoover baths). This is confirmed by the plots in the figures: the plot for 7 = 2 
reaches the asymptotic limit fastest, whereas 7 = 0.4, 10 have apparent exponents 
closer to 0.4 for small N. 

Thus the simulations of Mai et al. seem to give good evidence for a = 1/3 in the 
quartic FPU model and hence gives support for the ide a of a sin gle universality 
class. A discussion on these results is contained in refns. [H^, Ea). We note that 
the new prediction of a = 1/2 from MCT appears to be even harder to verify from 
simulations. 

Alternate mass FPU chain: Further support for the value a = 1/3 in the 
FPU system an d its universality comes from recent simulations of an alternate 
mass FPU chain [l57j . In this model one considers a chain with masses of particles 
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Figure 13. Plot of the heat current J versus system size for the alternate mass FPU chain, for different 
values of v = and with the mass ratio A = 1.5 (a). F ig (b ) shows the same data plotted with a scaled 
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at all even numbered sites being m\ and at odd numbere d site s being 1712 with their 
ratio being A = mi /mj . This system was first studied in 1181 ] where it was noticed 
that the temperature profile showed peculiar oscillations whose amplitude seemed 
to decay as N^ 1 / 2 with system size and scale linearly with the applied tempera- 
ture difference. This can be seen in Fig. (jlip . At the hotter end the lighter particles 
have a higher kinetic temperature, w hile at the colder end, the heavier particles 
are hotter. It was pointed out in refn. |l!8l | that the temperature oscillations make 
it difficult to define a local equilibrium temperature even at a coarse grained level. 
Temperature oscillations can in fact be seen even in an ordered binary mass har- 
monic chain but there seem to be some significant differences. In the harmonic 
case, there is a big difference between the case where iV is even and that where 
N is odd. This can be seen in Fig. (jlip where we plot the temperature profiles for 
chains of length iV = 128 and N = 129 for both the FPU chain and the harmonic 
chain with m\ = 0.8 and m<i = 1.2. For the harmonic case the oscillations for even 
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N are large and do not decrease with system size while those for odd N decrease 
with system size. In the FPU case there is not much difference between a chain 
with odd or even number of particles. Also for the harmonic case, in the bulk, the 
heavier particles are always hotter. 

However, even though there appears to be a problem in defining a local tempera- 
ture, one can still measure t he sy stem size dependence of the current in this system 
and this was done in refn. 157J ]. They studied the size dependence of the current 



for different values of the mass ratio A, keeping the average mass (mi + 772-2) /2 
constant. Remarkably it was found that at large enough system sizes the currents 
for different A all tend to converge to the same value. This can be seen in Fig. (|12p 
where one can see that the exponent is again as that at A = 1, i.e. a « 0.33. In 
this paper the authors next took a fixed value of A = 1.5 and studied the effect of 
changing the interparticle interaction strength (denoted as v = k^ in the paper). 
These results for current for different system sizes are shown in Fig. (|13j) . For small 
system sizes, one sees a flat region which is expected since for system sizes much 
smaller than the phonon-phonon scattering length scale, the system will behave 
like a harmonic chain. The scattering length should be larger for smaller v and 
this can be seen in the plot. At larger system sizes, all the curves tend to show the 
same decay coefficient with a ~ 0.33. A nice collapse of the data was obtained by 
scaling the system size by a length factor l{y) and this is shown in Fig. (113b ). The 
independence of the length parameter seems to be given by £(u) = l/tanh(2i/). A 
surprising point is that for any fixed system size, the value of the current saturates 
to a constant non-zero value as v — > 00. 

Double-well potential: It is interesting to consider simulation results obtained 
for the FPU interaction potential U(x) = k2X 2 /2-\-k^x /4 with neg ative k?. in which 
case we have a double- well potential. This case was first studied in 1351 ] which had 
initially reported a finite conductivity for this model but later it was found to 
have a power law divergence 0]. Here we present some new simulation results for 
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this model which in fact show that this model exhibits a fast convergence to the 
asymptotic regime with exponent a = 1/3. The parameter values hi = — 1, fej = 1 
and Tl = 1.25, Tr = 0.75 were considered. The temperature profiles in this system 
for different system sizes are shown in Fig. (114ft and are similar to the FPU profiles 
[Fig. ([!])] except that the boundary jumps are smaller. The inset of the figure shows 
a plot of J as a function of N . One can see that that by around iV = 512 the curve 
has reached the expected asymptotic slope corresponding to a = 1/3. Thus we 
again see evidence for a = 1/3. 

Discussion: In the absence of a rigorous proof it is fair to say that the question 
of universality of the heat conduction exponent and its precise value, in momentum 
conserving interacting systems in one dimensions, is still an open problem. This is 
especially more so since all the analytic methods and the exact result discussed in 
Sec. (|4.ip rely on use of the Green-Kubo formula for a closed system. As pointed out 
in sec (|2.2p the use of this formula and interpretation in systems with anomalous 
transport is not clear. The simulation results that we have presented in this review 
strongly suggest a single universality class, with a = 1/3, for momentum conserving 
interacting systems in ID. However one should probably bear in mind that in 
simulations, one can never really be sure that the asymptotic system size limit has 
been reached. It is possible for exponents to change in unexpected ways when one 
goes to larger system sizes. 

4-2.2. Momentum non- conserving models 

We will now look at heat conduction in one dimensional chains where the par- 
ticles experience, in addition to interparticle interactions, also external potentials 
which physically can be thought of as arising from interactions with a substrate. 
One of the first verification of Fourier's law in comp uter s imul ations was obtained 



by Casati et al. in the so-called ding-a-ling model lOll . Il22j ]. In this model one 
considers a system of equal mass hard point particles which interact through elastic 
collisions and where alternate particles are pinned by harmonic springs placed at 
fixed distances. The particles in between the pinned ones move freely. Clearly mo- 
mentum is not conserved and the authors, by studying system sizes upto N = 20, 
found evidence for diffusive behaviour. They calculated the thermal conductivity 
using both nonequilibrium simulations as well as using the Kubo formula and found 
good agreement between the two. We note that the system sizes studied in this 
paper are clearly too small to arrive at definite conclusions. Larger system sizes 
with the same parameter values were studied later by Mimnagh and Ballentine 
123]. They found that in fact the conductivity again started to increase as one 



went to larger sizes. Finally though, at system sizes N ~ 400, the conduc tivity 



again saturated at a new value which is much higher than that obtained in 1221 ] . 
This example nicely demonstrates the need for caution in drawing conclusions from 
small size data (also see discussion in 0] on these results). 

Since the work of Casati et al. , a number of papers have looked at heat conduc- 
tion in various momentum non-conserving models in one dimension and have all 
found evidence for the validity of Fourier's law. A model similar to the ding-a-ling 
is the ding-dong model and has all particles connected to fixed harmonic sp rings . 



This was studied by Prosen and Robnik and also shows Fourier behaviour 1241 ] . 
One of the first papers to recognize the fact that momentum non-conservation is 
a necessary condi tion to get finite heat conductivity in one-dimensional systems is 
that of Hu et al. |l25l |. From their simulations with various forms of Hamiltonians 
including, both a harmonic interparticle potential U (x) and a periodic onsite poten- 
tial of the Frenkel-Kontorva form ]V(x) ~ cos(ax)], they found that the presence 
of an external potential ty picall y led to a finite conductivity. The Frenkel-Kontorva 
model was also studied in |l2d ] who arrived at similar conclusions. A study of the 
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Fig ure 15. The geometry of the momentum-nonconserving, alternate-mass hardcore model, studied in 
[134|1 . The elementary cell (indicated by two dotted lines) has unit length 1 = 1. The bars have mass mi, 
and the particles have mass m.2- The two heat baths at temperatures and Tr are indicated. 



<f> 4 model [where U(x) = Jt2X 2 /2, V(x) = Xx 4 /4] by Hu et al. 1281 ]. again led to the 
conclusion of a finite conductivity. This study also emphasized the following point. 
Nonlinear integrable models usually give a flat temperature profile and making 
them non-integrable leads to a temperature gradient. However this in itself is not 
a sufficient condition to give a finite conductivity. 

Another nice simulation demonstrating the role of an external potential in giving 
rise to Fourier behaviour is that on the binary mass hard particle gas [HI. Mo- 



mentum non-conservation is ensured by confining the motion of alternate particles 
inside finite cells while allowing them to interact, through elastic collisions with 
neighbors [shown in Fig. (|15j) ]. The nonequilibrium simulations in this paper with 
Maxwell heat baths (with N ~ 512) convincingly shows the validity of Fourier's 
law and also the presence of local thermal equilibrium. Secondly, equilibrium simu- 
lations were also performed to compute current-current correlation functions, and, 
using the Green-Kubo formula a value of k close to the nonequilibrium result was 
obtained. It is worth noting that this model has zero Lyapunov exponent and thus 
is non-chaotic. A related study is that in refn. 1331 ] who studied heat conduction 
in a gas of hard rods placed in a periodic potential. 

Results for the <j) 4 model: We will describe in some more details work on the 
(ft 4 model [U(x) = k2X 2 /2, V(x) = Xx 4 /4\ which appears to be one of the most 
well-studied of the momentum non-conserving models and where some analytic 
results have also been ob tain e d. H eat conduction in the (ft 4 model was first studied 
by Aoki and Kusnezov 1271 . Il29f ] who performed both nonequilibrium measure- 
ments as well as Green-Kubo based equilibrium measurements. Studying system 
sizes upto N = 8000 they concluded that this system had a finite conductivity 
and Fourier's law was valid. The value of k, obtained from the nonequilibrium 
measurements and from the Green-Kubo formula were again shown to be in good 
agreement. The authors also numerically obtained the temperature dependence of 
k and found k(T) ~ j 1-1 - 35 . A number of other papers have performed simulations 
of the (ft 4 mode l and studied various aspects such as the spreading of localize d dis - 
turbances [l28l ] and the dependence of thermal conductivity on temp erature [l32 |. 
The model was studied analytically by Lefevere and Schenkel |l3d( ] and later by 
Aoki et al. |l3ll ] using a Peierls-Boltzmann kind of approach for the case of weak 
anharmonicity and they too obtained a finite conductivity. They however obtained 
a temperature dependence k ~ 1/T 2 and this is probably the correct low tem- 
perature (corresponding to weak anharmonicity) behaviour, since kinetic theo ry is 



expected to be reliable in this regime. Direct nonequilibrium simulations in [131 ] 
infact found reasonable agre emen t with the predictions from kinetic theory, at low 
temperatures. The study in 1321 ] however finds a somewhat different temperature 
dependence at low temperatures (k ~ 1/T 1 ' 56 ). We note that a scaling property of 
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Figure 16. Plot of the temperature profile in the </> 4 model for different system sizes. The inset shows the 
system size dependence of the current and the straight line drawn has slope —1. Error bars are much 
smaller than symbol sizes. The parameter values for the figure are fc2 = 1, A = 1,Tl = 1.25, Tr = 0.75 

and bath coupling 7 = 1. 



k(T, A), to be discussed later in Sec. ([5]), implies that k = k(XT). 

In Fig. (|16p we show typical plots of the temperature profile in the <^> 4 chain. The 
inset in the figure shows the 1/N dependence for the current. These simulations 
were performed using white noise Langevin dynamics using the velocity- Verlet 
algorithm and for the largest system size (N = 8192) required ~ 2 x 10 9 time steps 
with At = 0.0025, to equilibrate. 

Before concluding this sectio n we mention the results on a clas s of rotor mod- 
els studied by Giardina et al. 135| and Gendelmann et al. 136 ]. These models 



were originally proposed as examples of momentum conserving systems which gave 
Fourier behaviour. The interparticle potential is taken to be U{x) = 1 — cos(x) and 
the onsite ter m V(x) = 0. Nonequilibrium and equilibrium (Green-Kubo based) 
simulations in Il35 |. Il36j ] both indicated that this model gave a finite conductivity. 
The paper by [1361 ] also reports a phase transition from infinite to finite conduc- 
tivity as a function of temperature. Given that these simulations are restricted to 
relatively small sizes (upto N = 2400), one su spects that thi s is probably a cross- 
over effect. Simulations for larger systems in 137 . 1381 . 158 ] indeed suggest that 
there may not be any true transitions and that, at all temperatures the asymptotic 
beaviour is Fourier-like. An analytic study of the rotor model using self-consistent 
reservoirs ( wit h vanishingly small coupling to interior points ) has also claimed a 
transition |l39l | . A similar claim of possible transitions from finite to diverging con- 
ductivity in other momentum no n-con serving models such as the Frenkel-Kontorva 
and 4> 4 model has been made in [l4d( |. 

The fact that a momentum conserving model gives finite conductivity is at first 
surprising. However given the form of the interparticle potential in the rotor model 
it is probably more physical to think of this model as an angular momentum 
conserving model rather than linear momentum conserving one. Thus it seems 
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more natural to think of the position variables xi as transverse angular degrees 
of freedo m. I n this case one expects different hydrodynamic equations (see for 
example 14l| ) and the Fourier behaviour observed is then not surprising. One can 



also think of the rotor model as the classical limit (large spin) of quantum lattice 
spin chain models which also are momentum-nonconserving. 

Quantum mechanical models: The study of nonequilibrium steady states of 
interacting quantum systems by simulations is an important and difficult problem. 
There have been a few attempts at addressing th is issue, a nd we will summarize 
these. The first set of papers were by Saito et al. 1431 . 144], who used the master 



equation approach to connect different temperature reservoirs to a quantum spin 
chain. One interesting result was that a temperature gradient was formed for the 
case where the model was non-integrable, while a flat profile was obtained for an 
integrable model. A study of the current-curent corr elator yielded a power law 



decay C{i) ~ f~ 1,5 implying a finite conductivity 14411. 

In another interesting work, Mejfa-Monasterio et al. [l45| | have devised what they 
call a quantum stochastic reservoir. Using this they have performed nonequilib- 
rium simulations, again of a quantum spin chain. They also observe a temperature 
gradient for the non-integrable model and a flat profile for the integrable model. 
Further they measured the nonequilibrium steady state current for different sytem 
sizes and found a J~iV _1 dependence for the non-integrable case and J ~ N° 
for the integrable case. 



5. Systems with disorder and interactions 

As discussed in sec. (j3.4|) localization of eigenfunctions or of normal modes strongly 
affects transport in materials containing random impurities. In electronic sys- 
tems localization has its strongest effect in one dimensions where any finite dis- 
order makes all eigenstates localized and one has an insulator. The presence of 
inelastic scattering, such as is caused by electron-phonon interactions, leads to 
hopping of electrons between localized states and gives rise to a finite conduc- 
tivity. The question as to whether electron-electron interactions lead to a simi- 
lar effect has attracted much att ention recently but is still not fully understood 

m, H EH, EH, EH, EHU, H . The main interest is to understand the transi- 



tion, from an insulating state governed by the physics of Anderson localization, to 
a conducting state as one increases interactions. One can ask the same question in 
the context of heat conduction by phonons and consider the effect that phonon- 
phonon interactions have on localization. Here we will mainly discuss the effect of 
anharmonicities on the steady state transport of heat through a chain of oscillators 
with random masses . The effect of interactions between phonons on localization 
caused by di s order has also been investigated by looking at the spreading of wave 



packets 15ll . Il52l . Il59l . Il60l ] and we will briefly discuss these results at the end of 
this section. 

An early work on steady state heat cond uction in disordered anharmonic systems 
is that of Payton, Rich and Visscher |l53j ] who studied mass-disordered lattices in 
the presence of cubic and quartic interparticle anharmonicities. They performed 
nonequilibrium simulations with stochastic baths in one and two dimensions. Their 
main conclusion was that in most cases interactions (interparticle anharmonicity) 
seemed to greatly enhance the conductivity of the system (except for the case of 
very weak disorder). We note that, at that time simulations were restricted to small 
sizes and it was wrongly assumed by the authors that the disordered harmonic 
lattices in one and two dimensions, as well as the anharmonic ones, had finite 
thermal conductivities. The system size dependence was not studied systematically. 
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Similarly a study by Poetzsch and Bottger for a two dimensional lattice 

system found that, while quartic anharmonicity enhances the conductivity of the 
disordered system, cubic anharmonicity reduces it. Again this study was restricted 
to small system sizes and assumed that the conductivity is finite. 

The first systematic study of the joint effects of anharm oncit y and disorder on 
the system-size dependence of heat current was by Li et al. 1551 ]. They studied the 



mass-disordered FPU chain using Nose-Hoover nonequilibrium simulations. Their 
conclusion was that this model showed a transition, from a Fourier like scaling 
J ~ N^ 1 at low temperatures, to a pure FPU like behaviour with J ~ jV~ ' 57 
at hig h te mperatures. A more recent simulation of the same model by Dhar and 
Saito [l57 ] suggests that this conclusion may be incorrect. What Li et al. observe 



is probably a cross-over effect and there is really no true transition in transport 
properties. It appears that disorder becomes irrelevant as far as the value of the 
exponent a is concerned. We will here present the latest simulation results of the 
disordered FPU chain as well as result s from the study of the disordered <p 4 lattice, 
where similar conclusions are reached [1561 ] . Temperature driven phase transitions 
in one dimensional heat conduction have also been reported in some other models. 



However as has been pointed out in refn. [1581 ] these are probably cross-over effects 



and there is no true transition in these models. 

The general form of the Hamiltonian that has been studied by various people is, 
in the one-dimensional case, given by: 



9 2 4 

^ v 2mi 2 4 J 

l=l,N 

+ £ [ k {Xl ~ X 2 l ~ 1? + v {Xl ~ ] (118) 

l=l,N+l 



with fixed boundary conditions xq = xn+i = 0. The masses {m{\ are chosen 
independently from some distribution p(m), e.g. one uniform in the interval (m — 
A, rh + A) or a binary distribution given by P{m) = 5[m — {in — A)]/2 + 5[m — 
(m + A)]/2. The chain is connected at its ends to two heat baths at temperatures 
Tl and Tr respectively. Here we will mostly consider white noise reservoirs, but 
will also give some results with Nose-Hoover baths. The equations of motion of the 
chain are then given by: 

mixi = -k Q xi - Ixf - k(2xi - - xi + ±) 

- u[(xi - x ; _i) 3 + (xi - x t+1 f] - -yi±i + 7]i , (119) 



with rji = r/^e^i + rmS^N, 7; = 7(^,1 + <5«,tv)i an d where the Gaussian noise 
terms satisfy the fluctuation dissipation relations {r)L(t)r]L(t')) = 2 r ykBTi,5(t — t'), 
{r] R (t)VR(t')) = ^k B T R 5(t-t'). 

Note that Eq. (jll9H is invariant under the transformation Tl r — ► sTl^r, {x{\ — > 
{s l l 2 xi} and (I, v) — > (l,u)/s. This implies the scaling relation J(sTl, sTr,1,v)) = 
sJ(Tl,Tr, si, sv). For the conductivity k this implies n = k(vT, IT). Thus the effect 
of changing temperatures can be equivalently studied by changing anharmonicity. 
We will first discuss the unpinned (momentum conserving) case and then the pinned 
(momentum non-conserving) case. 

Disordered FPU chain: This corresponds t o ta king k Q = A = in the Hamil- 
tonian in Eq. fjl 18H . and is the case studied by |l55l ] and in 157 ]. There are two 
important parameters in the problem, namely the disorder strength given by A and 
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Figure 17. Plot of heat current versus system size, for th e dis ordered FPU chain, for different values of 
v. The data in the inset corresponds to parameters as in [155| . namely (7l,T^) = (0.001, 0.0005) with 
Gaussian white noise bath for v = 1 (WN) and v = 0, and Nose-Hoover bath (NH) for v = 1 (from |157|| ). 



the anharmonicity given by v. Let us consider the two limiting cases, of the disor- 
dered harmonic chain (y = 0,A/ 0), and of the ordered FPU chain (y ^ 0, A = 0). 
For the former, with fixed boundary conditions it is expected that a = —1/2, while 
for the ordered FPU chain one expects a = 1/3. In the presence of both disorder 
and interactions a possible scenario is that for strong disorder one gets a = —1/2 
while with strong interactions, one gets a = 1/3 and there is a phase transition 
between the two behaviours as we change parameters. The numerical results that 
we will discuss, suggest that there is no such transition. Note that in both the 
limiting cases, the low frequency long wavelength modes are believed to play an 
important role in tra nspo rt. 

The simulations in [1571 ] looked at the case of binary mass distribution with fh = 
1, A = 0.2 and different values of the interaction strength v = 0.004, 0.02, 0.1, 2.0. 
Averages were taken over 50 — 100 samples for N < 1024, 10 samples for N = 
1024 - 16384, and 2 samples for N = 32768 and 65536. In Fig. (JT7|) the results of 
simulations for the disorder averaged current [J] for v = 0.004, 0.02 and v = 0.0 
are shown. For small values of v one sees that, at small system sizes the current 
value is close to the v = value. As expected one has to go to large system sizes to 
see the effect of the weak anharmonicity. At sufficiently large the same system 
size dependence of J is obtai ned as that for the ordered FPU chain, namely with 
a = 1/3. The authors in [157] then show that by scaling the current by appropriate 
factors, the data for the disordered case can be made to collapse on to the binary- 
mass ordered case. This is shown in Fig. (]18p (for v = 0.02,0.1,2.0). Thus these 
results show that the asymptotic power law dependence of the current is always 
dominated by anharmonicity and there seems to be no transition. Disorder only 
decreases the over all c onductance of a sample. 

The authors of |l57| have also investigated the parameter range studi ed in 1551 ] 
and explained the reasons which led to the erroneous conclusions in [l55 |. of a 
transition in conducting properties at low temperatures (or equivalently small an- 
harmonicity). In fact this can be understood even from the data for [J]N in Fig. (I17p 
for v = 0.004. We see that at around N ~ 1000 — 2000 the data seem s to flatten 
and if one had just looked at data in this range, as was done by [l55l ]. one would 
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Figure 18. Plot of scaled heat current [J] s for the disordered FPU chain and the curre nt J for the 
ordered chain, for different values of v. The x— axis is scaled as in Fig. H13H (from [157(1 ^ . 



conclude that Fourier's law is valid. However the behaviour changes drastically 
when one looks at larger system sizes and one again gets the usual FPU behaviour. 
The inset of Fig. (|17p shows results for parameters as in 155] but for much larger 
system sizes. This case corresponds to a much smaller value of v and so it is ex- 
pected that it will follow the v = curve till very long length scales and this is 
clearly seen. However, at around N = 16384, there is a tendency for the curve to 
turn up and it can be expected that the same asymptotic behaviour to eventually 
show up. While a transition cannot be ruled out at even lower temperatures and 
with stronger disorder, this seems unlikely. Also, if there is such a transition, it 
should probably be to a disordered phase with [J] ~ iV~ 3 / 2 . 

It is interesting to consider the temperature dependence of conductivity in the 
disordered FPU chain. The scaling property of the current, mentioned earlier [after 
Eq. (|119p ]. implies that the thermal conductivity has the form n = k(vT). For small 
anharmonicity (y « 1), the earlier results for the ordered alternate mass FPU 
chain imply that at large system sizes k ~ iV 1 / 3 jv 2 ^ and from the scaling property 
this immediately gives k ~ 1/T 2 / 3 at low temperatures. However at small system 
sizes [N << we expect the system to behave like a harmonic system with 

k ~ T°. At high temperatures the conductivity will saturate to a constant value. 
Experimentally, the temperature dependence of the thermal conductivity may be 
easier to measure and one can verify if this is unaffected by disorder [see, for 
example sec. (JSJ)]. 

Disordered 4 chain: Let us now look at the case where the particles are 
subjected to an external pinning potential in addition to nearest neighbor harmonic 
interactions. We will consider the anharmonicity to be an onsite quartic term (thus 
A > 0, v = 0, also k , k > ) in which case this corresponds to the discrete 4> 4 model. 
Pinning greatly enhances the difference between heat transport in a random chain 
with and without anharmonicity and thus is a good testing ground for the effect of 
anharmonicity on localization. This model is also closer in spirit to charge transport 
by hopping in random media. Again let us look at the two limiting cases. In the 
case with a pinning potential at all sites the disordered case (A = 0,A / 0) gives 
J ~ e~ cN . For A = and A / we have seen from sec. (|4.2.2p that one expects 
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Fourier's law to be valid and so a = 0. 

The case with parameters k = k Q = 1 , A > and a uniform mass distribution 
with fh = 1.0 and A = 0.2 was studied in [il56;]. In Fig. (|19h the result of simulations 



for different values of anharmonicity A = 0.004 — 1.0 is given. As can be seen from 
the data, there is a dramatic increase in the heat current on introduction of a small 
amount of anharmonicity and the system-size dependence goes from exponential 
decay to a 1/N dependence implying diffusive transport. For smaller A the diffusive 
regime sets in at larger length scales but as in the FPU case, here too one finds that 
anharmonicity determines the system size scaling and no transition is observed. A 
measure of the relative strengths of anharmonicity and disorder is obtained by 
looking at the ratio of the energy scales E a = l(x )/4 and E& = TA/m. For the 
given parameters one finds e = Ea/Ej ~ 0.008 for A = 0.004. Unlike the FPU case, 
in this model, it does not seem that any simple scaling of the data is possible. 

Thus this study shows that introduction of a small amount of phonon-phonon 
interactions in the disordered harmonic chain leads to diffusive energy transfer, i.e., 
the insulating chain becomes a normal heat conductor. How exactly this occurs is 
not clear. It is possible that anharmonicity gives rise to extended states or leads 
to hopping of energy between states which are now approximately localized (i.e 
they are no longer exact normal modes, but have a small rate of energy leakage 
to nearby modes) . There is no evidence of the existence of a finite critical value 
of anharmonicity required for this transition. For small values of anharmonicity it 
is necessary to go to larger system sizes to see the transition from insulating to 
diffusive. As in the FPU transition to a localized phase at a very small 

value of anharmonicity is possible and would be difficult to observe in simulations, 
because equilibration times increase rapidly with decreasing A. 

An interesting question in this model is the limiting behavior of k(XT) for 
(AT) — > 0. It turns out that the temperature profiles for the disordered </> 4 chain are 
qualitatively different from the ordered chain and this means that the temperature 
dependence of conductivity is different for the two ca ses. For the ordered case, from 
kinetic theory one gets k ~ 1/(A 2 T 2 ) for small AT [1311 ] . while for the disordered 
case |l56l ] found k ~ (AT) 1 / 2 . For the FPU chain on o ther hand, the ordered and 
disordered cases give similar temperature profiles [l57j |. 

There have been a number of studies on disordered anharmonic chains which 



ected into a system 
looked at FPU 



have investigated the spreading of localized pulses of energy in 
at zero temperature. The study by Bourbonnais and Maynard 
type of systems in one and two dimensions and observed that anharmonicity desta- 
bilizes the localized modes and the diffusion of pulses was found to be anomalous. 
This seems to be consistent with the heat conduction results on the disordered 
FPU chain. A similar zero temperature stud y of the mass disordered FPU system 
was carried out by Snyder and Kirkpatrick [160f | who however found evidence for 
diffusive transport at sufficiently strong anharmonicity. In the case of the </> and 
related models there have been some extensive recent numerical studies and here 
the conclusions are somewhat contradictory to the heat conduction results. The 



sprea ding of localized energy pulses has been reported to be su b-diffusive in [151 ] 



while 1521 ] reports absence of diffusion. The study in refn. 152 ] offers a picture of 



spreading of an initially localized energy wavepacket to a limiting profile as taking 
place through nonlinearity induced coupling between the localized modes. All these 
studies suggests that the behaviour of a heat pulse at zero temperature a nd t hat 
at finite temperature can be very different. Indeed as pointed out nicely in [1421 ] . it 
is necessary to look at appropriate spatiotemporal correlation functions of closed 
systems at finite temperatures in order to understand diffusion in the open system. 
This of course is also what one effectively does in the Green-Kubo approach. 
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Figure 19. Plot of [J]N versus TV for the disordered (f> 4 chain, for different values of A. The inset show s 
results obtained for the case with interparticle anharmonicity and onsite harmonic pinning (from 156]). 



Finally we note an interesting related problem that was studied by Rich and 
Visscher [65] . They considered heat conduction in a disordered Harmonic chain with 
self-consistent reservoirs. Since self-consistent reservoirs can be roughly though of as 
some sort of nonlinear ity leading to incoherent scattering of phonons, this problem 
has some similarity with that considered in this section. Based on exact numerical 
calculations on small chains, their main conclusion was that the presence of self- 
consistent reservoirs leads to a finite conductivity for chains with both free and 
fixed boundary conditions (and no bulk pinning). The self-consistent reservoirs 
makes the model momentum non-conserving so this is consistent with the results 
of the disordered 4 chain presented here. A very interesting conjecture made in 
this paper is that a finite conductivity will be obtained if the limits N — > oo 
first , and then coupling to self-consistent reservoirs — > are taken. A recent paper 



1611 ] has studied heat conduction in a disordered harmonic chain with an energy 



conserving stochastic dynamics and has obtained rigorous results which indicate a 
finite thermal conductivity of the system. 



6. Interacting systems in two dimensions 

We have seen that, in the one-dimensional case, it is usually quite difficult to obtain 
the asymptotic system size dependence of the current. In order to get the correct 
exponent requires one to go to large system sizes and at some point the sizes 
required are beyond current computational capabilities. Of course a combination 
of simulations and results from analytic work gives one some confidence about 
the results obtained so far. In the case of higher dimensional systems naturally 
one can expect the same computational difficulties and in fact here they become 
more pronounced since the number of particles is now L d where L is the linear 
size and d the dimensionality. The good news is that there is general agreement 
on the system-size dependence of conductivity from different analytical methods. 
Both MCT [3] and the hydrodynamics approach [72] predict that for a momentum 
conserving system, the thermal conductivity diverges logarithmically with system 
size in 2D and is finite in 3-D. In the presence of pinning all theories predict a finite 
conductivity in all dimensions. 
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There are few simulation studies in higher dimensional systems and we summa- 
rize the main results obtained so far. Early studies of heat conduction were mainly 
interested in finding the te mper ature dependence of thermal conductivity and as- 
sumed that this was finite [\M. 
size dependence was probably that by Jackson and Mistriotis [163]. They studied 
the diatomic Toda lattice and concluded that the thermal conductivity was finite 
for mass ratio greater than a critical value and diverged otherwise. 

Mo re e xtensive studies on the system size dependence were made by Lippi and 
Livi Hfii| for an oscillator system in two dimensions with vector displacements 
and interparticle interactions. They considered a L x x L y lattice with the following 
Hamiltonian: 

ff = EEil +f/ (i x '+i.i- x .ii)+%,i+i-si) , (120) 

i=l j=l 



162, 1641. One of the first paper to study system 



where Xjj denotes the vector displacement (taken to be two-dimensional vectors) 
of a particle at lattice site where i = 1,2..L X and j = l,2...L y and py de- 
notes the corresponding momentum vector. Two kinds of interparticle potentials 
were studied, namely a FPU type potential given by U(x) = x 2 /2 + k^x^/A and a 
Lennard-Jones potential given by U(x) = A/x 12 — B/x 6 . Both models gave similar 
results. Nonequilibrium simulations using Nose-Hoover baths, as well as equilib- 
rium simulations based on the Kubo formula, were performed. Nonequilibrium 
simulations were first performed on strips of width L y with aspect ratio L y /L x < 1 
and with heat conduction in the x-direction. It was observed that for fixed L x , 
as one increased L y , the current seemed to saturate to a constant value for quite 
small values of L y /L x . Subsequently, to save on computational time, the authors 
considered the value L y /L x = 1/2 in all their simulations. Studying system sizes 
upto L x = 128 they obtained a logarithmic divergence, with system size, of the 
conductivity i.e. k ~ ln(L x ). The equilibrium simulations, performed over similar 
system sizes, and using a microcanonical ensemble gave a t~ l dependence for the 
current-current correlation function. Using the Green-Kubo formula this implies 
again a logarithmic divergence, with system size, of the conductivity. 

The same Hamiltonian as in Eq. (|120p with FPU interactions but with a scalar 
displacement field was studied later by Yang and Grassberger |l66l ]. This paper 
looked at somewhat bigger system sizes than 1651 ] but were unable to verify the 
logarithmic divergence and instead obtained a power law dependence with an ex- 
ponent a ~ 0.22. A careful investigation of the value of r = L x /L y , at which a 
dimensional cross-over from ID to 2D behaviour occured was carried out. Their 
conclusion was that at large values of r, the conductivity k diverged as a power 
law with a = 0.37 ± 0.01 while for small r they obtained a ~ 0.2. The data for 
conductivity versus system size for different values of r is shown in Fig. (|20p . The 
conductivity plotted in the figure was defined as 



^center — , irp i , s > (l^l) 

yul j dX ) center 



where (dT/dx) cen ter is the temperature gradient evaluated numerically at the cen- 
ter. This definition was used to take care of the boundary temperature jumps that 
are usually present for small system sizes [see Fig. (|21h ]. Based on the data in 
Fig. (|20p the authors also made the interesting suggestion that the cross-over from 
ID to 2D behaviour takes place at r — > oo in the limit L x — > oo. This has obvious 
implications for experimental tests on the dependence of conductivity on length, 
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Figure 21. Temperature profiles for scalar 2D FPU lattices with L x = 
temperatures of the heat baths at both ends were fixed at (Tj,, Tjj) 



128 and L y = 1, 3, and 32. The 
= (10.0, 6.0) (data from [161 ). 



for systems such as nanotubes and nanowires. Another paper [l67l ] again study- 
ing the vector model for even larger system sizes (upto 64 x 65536) has claimed 
observing a logarithmic divergence. However one of the authors of the pape r has 
expressed doubts about whether convergence has been attained at these sizes |l68j | 
and this seems very likely to be the case. 

The most recent simulations by Shiba and Ito 170\ consider ed th e same Hamil- 
tonian as in Eq. (|120j) and used the same parameter set as 1651 ]. namely k± = 
0.1, Tl = 20, Tr = 10. They performed nonequilibrium Nose-Hoover simulations 
and studied system sizes upto 384 x 768. Their data for conductivity versus system 
size is plotted in Fig. (1221) . The exponent a ~ 0.268 obtained by them appears to 
be significantly different from logarithmic behaviour. We also show the tempera- 
ture profiles for different system sizes [Fig. (|23j) ] and it appears that the boundary 
jumps are quite negligible. 
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Figure 22. System size dependence of the thermal conductivity for 2D FPU lattice with L y : L x = 1:2, 
plotted on a log — log scale. The dashed line represents the result of a power-law fitting in the region 
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Figure 23. Temperature profiles for the 2D FPU lattice with L y : L x = 1:2. The sequences represent 
the results for the sizes L x = 192, 384, and 768. The temperatures of the heat baths at both ends were 
fixed to (Tj,,Tfj) = (20.0, 10.0). The horizontal axis represents the position in the x-direction, scaled by 
the system size L x , and the vertical axis represents the local temperature (data from [170j | V 



Unlike in ID, where the hard particle gas has been intensely studied, there have 
bee n ver y few studies on hard disc systems. For a hard disk fluid system, Shimada et 



al. [IflU reported a to be less than 0.2. Thus for momentum-conserving systems in 



2D it is fair to say that simulations have not been able to unambiguously establish 
the logarithmic divergence of the conductivity predicted from theory. Further work 
is clearly needed here. 

As far as momentum non-conserving interacting systems are concerned one would 
naturally expect Fourier behaviour, given that this is the case even in one dimen- 
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sion. In the next section we will discuss a number of momentum non-conserving 
models of non-interacting particle systems which can be shown (including some 
rigorously) to satisfy Fourier's law. In these models noninteracting particles are 
scattered from fixed scatterers. These models however suffer from the drawback 
that there is no mechanism for local thermal equilibration and so the meaning of 
temperature and Fourier's law is somewhat artificial. There have been a few papers 
which have introduced particle interactions in these kind (hard particle scattering) 
of models. Here we discuss two such models. 
The first model, introduced by Mejia-Monasterio et al. [l71 . 1721 ] . is one in 



which noninteracting particles move among a periodic array of circular scatterers 
[see Fig. (|24p ]. The dynamics is specified as follows. Consider dimensionless units 
such that the mass of the moving particles is 1 and the moment of inertia of 
the scatterers is rj. Then, if (v n ,Vt) are the normal and tangential components of 
velocity of the particle at the time of collision, and to is the angular velocity of the 
discs, then after the collision they are transformed to (y' n ,v' t ,to') which are given 
by the linear transformation: 



v'„ = -v. 



v t = v t ~ YT^j^ Vt ~ 

to' = to + ^— (v t -to) . (122) 

The dynamics conserves total energy + + nto 2 /2 and angular momentum 
and is time-reversal invariant (however, the transformation is non-syplectic). This 
system was then connected to two reservoirs of both heat and particle and which 
are specified by temperature and chemical potentials (Tl, /jll) and (Tr, /Ur) respec- 
tively. Thus both heat and particle currents were generated. Performing detailed 
simulations on this system, some of the main conclusions of the paper were: (i) the 
system satisfied local thermal equilibrium, (ii) both heat and particle currents sat- 
isfied usual linear response relations with finite transport coefficients, (iii) Onsager 
reciprocity relations were satisfied. The largest system studied had about 100 discs 
in the conducting direction (and two discs in the vertical direction). Note that in 
this model interactions between particles is introduced indirectly. Motivated by this 
model, refn. 173j ] studied an idealized model with noninteracting tracer particles 



moving between fixed energy storing centres and exchanging energy with these. Lo- 
cal thermal equilibration and temperature profiles were analytically studied in this 
work. Another model where an explicit verification, of linear response relations for 
heat and particle transport were obtained, as well as Onsager reciprocity relations, 
is a ID electronic system with self-consistent rese rvoirs |70l| . 
Another recent study by Gaspard and Gilbert [l74l . 1 1751 . 176j ] has considered a 



system where hard disc particles are confined within periodic array of cells formed 
by fixed scatterers. The model is explained in Fig. (|25|) , The main idea of the 
authors has been to introduce a three time-scale mechanism in generating the heat 
conduction state: (i) a short time scale r wa u corresponding to particles motion 
within a cell with negligible energy transfers, (ii) an intermediate time scale Tjyi nar y 
corresponding to binary collisions which lead to local equilibrium and (iii) a long 
time scale T macro of the macroscopic relaxation of Fourier modes. Based on a master 
equation approach the authors are able to demonstrate the validity of Fourier's law 
and find an explicit expression for the thermal conductivity of the system. Note 
the similarity of the model with the one dimensional model described in Fig. (I15h . 
though the mechanism leading to Fourier behaviour there is possibly different. 
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Figure 25. Model studied in [17411 . The small colored discs move among a periodic array of fixed black 
discs. Each small disc is confined to move in a cell bounded by four fixed discs. Most of the time the 
small disc moves with constant energy and undergoing elastic collisions with the fixed discs. Once in a 
while there is a collision between particles in two neighboring cells and there is exchange of energy. The 
solid lines show the trajectories of the centers of the moving particles about their respective cells. The 
colors are coded according to the particles kinetic temp eratu res (from blue to red with increasing 

temperature) (from [17410 , 



As far as simulations of oscillator systems is concerned, a t wo d i men sional 
momentum non-conserving lattice model was studied by Barik |l77 . 1781 ] who 
studied a scalar displacement model with harmonic interparticle interactions and 
an onsite potential V(x). Several different forms of V{x) were studied. With 
a Prenkel-Kontorva interaction given by V n ( x) = cos(x) the author reports a 



power law divergence of the conductivity 1771 ] . For two other models, of the form 
Vb(x) =_— cos(x) — sin(2x)/2, a logarithmic divergence of the conductivity was 
found [l78 |. However the reason why a finite conductivity has not been obtained 
in these momentum non-conserving systems is probably because the system sizes 
studied are too small (upto 240 x 240). This is evident from the quite large bound- 
ary temperat ure jump s that can be seen in the temperature profiles reported in 
these papers 177 . 178| . This means that the contact resistances are contributing 
significantly to the measured resistance. As we have discussed earlier in sec. (14.2.1]) 
this will give a higher apparent divergence of the conductivity than is actually the 
case. In this case looking at ^center , defined in Eq. (I12ip . may be a good idea. 
Another point is that speci fic fo rms of interaction strengths also might lead to a 
faster convergence. Thus in [l78l ] it can be seen that for the same applied temper- 
ature difference and same system size, the potential V a (x) gives a flat temperature 
profile while with Vb(x) one gets a significant gradient. This feature has also been 
observed earlier in sec. (14.2, lh where we saw that the random collision model and 
the ID double-well potential gave fast convergence to the asymptotic limit. 

It appears that the simulation results in two dimensions are quite inconclusive 
for momentum conserving systems as regards the system size dependence of con- 
ductivity. More extensive studies with larger system sizes and different forms of 
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interaction potentials are needed to confirm the theoretical predictions. For mo- 
mentum non-conserving lattice systems it is likely that simulations on larger system 
sizes will verify validity of Fourier's law. 



7. Non-interacting non-integrable systems 

Probably the first rigorous demonstration of Fourier-like dep ende nce of the cur- 
rent in a Hamiltonian system was by Lebowitz and Spohn 1791 ] in the Lorenz 



gas model. This model consists of a gas of classical point particles in dimension 
d > 2 which undergo elastic collisions with fixed randomly placed spherical scat- 
terers. The authors studied this system with stochastic boundary conditions on 
two bounding walls corresponding to temperatures Tl and Tr. They could prove 
rigorously that in the Boltzmann-Grad limit of a large number of small scatterers, 
one gets J ~ (Tl — Tr)/L, where L is the length of the box. 

A number of quasi-l-D Lorenz-gas-like systems have been numerica lly i nvesti- 
gated recently and have provided some more insight. Alonso et al. 1801 ] stud- 
ied a channel with an array of periodically placed semicircular scatterers. From 
their nonequilibrium simulations with Maxwell baths, they verified Fourier's law 
and showed that the nonlinear temperature profile could be understood by a 
simple model of diffusing particles with an energy dependent diffusion constant 
D(E)~E 1 / 2 . 

The particle dynamics of the Lorenz-gas model with convex scatterers is known 
to be chaotic. Li et al. [l8lj ] explored the question as to whether a positive Lya- 
punov exponent, i.e. , a chaotic dynamics, is a necessary condition in order to 
get Fourier's law. They considered a system of non-interacting particles moving 
in a quasi ID channel and suffering elastic collisions with fixed triangular shaped 
scatterers placed throughout the channel. The cases of a regular, and a random, 
array of scatterers was considered. It can be shown that in both these cases, the 
dynamics has zero Lyapunov exponent unlike the case with convex scatterers. The 
authors found from their nonequilibrium simulations that while the regular array 
gave a diverging thermal conductivity, the random array gave a finite conductivity. 
In both cases a temperature gradient was obtained, though they had very different 
forms, highly nonlinear in the periodic case and linear in the random case. Thus 
the random array of non-chaotic scatterers gave rise to normal Fourier heat trans- 
port in the channel. Correspondingly it was shown that the particle motion was 
superdiffusive for the periodic case and diffusive for the random case. Thus this 
simulation shows tha t cha os is not a necessary condition for diffusive transport. 



In other studies [H3, E13] it has been seen that, even with a periodic array of 



triangular scatterers, one gets Fourier transpor t wh enever the internal angl es are 
irrational multiples of it. Finally Li and Wang 184 ] and Denisov et al. [l85 ] have 



given analytic arguments to relate the diffusion exponent of the heat carriers to 
the conductivity exponent a. That chaos is not a necessary condition for normal 
transport was also seen earlier in sec. (jl]). There we saw that the interacting hard 
particle dimer gas (which is non-chaotic) gives normal transport in the presence of 
an external potential. 

The validity of Fourier's law in Lorenz-gas models basically arises due to the 
diffusive motion of the heat carriers. However the absence of interactions makes 
these models s ome what ill-behaved from the thermodynamic point of view. As 
pointed out in [l~8^ these models lack local thermal equilibrium and so the mean- 
ing of temperature and Fourier's law in these systems is somewhat different from 
that one usua lly has in nonequilibrium thermodynamics. The models studied in 



17ll . Il72j . [174 ], and discussed in sec. ([6]), are examples of similar momentum non- 
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conserving models of particles with collisional dynamics where however the intro- 
duction of interactions leads to local thermal equilibrium. 



8. Experiments 

The experimental measurement of thermal conductivity of a system is much more 
difficult than, for example, its electrical conductivity. One realizes this from the 
simple fact that it is easy to construct an ammeter to measure electrical current 
but that is not so in the case for heat current. Thus measurements of thermal 
conductivity require special methods and often the interpretation of experimental 
data themselves require involved theoretical modeling. This is perhaps one reason 
as to why there has, until recently, not been much experimental studies which have 
addressed the precise question of the system size dependence of thermal conduc- 
tivity and its expected divergence in low dimensional systems. The situation has 
changed recently with the advent of nanophysics. Understanding heat transfer in 
systems such as nanowires and nanotubes is not only a question of basic interest 
but of technological importance too. With amazing advances in nano-technology 
it has now become actually possible to measure the thermal conductivity of a nan- 
otube suspended between two thermal reservoirs. Here we briefly discuss some of 
the experiments on nanowires and nanotubes. We will try to explain the present 
understanding and also try to emphasize the relevance, of the knowledge that has 
been obtained from studies of simple models discussed in this review. In all the 
experiments that we will describe here the heat current is believed to be mainly 
due to phonons. 

The most common approach to understanding experimental data on heat con- 
duction is perhaps through simple kinetic theory picture which says that the con- 
ductivity in a phonon system is proportional to cv£ where c is the specific heat per 
unit volume, v the sound speed and I the phonon mean free path. For system sizes 
smaller than £, one expects ballistic transport and roughly one can replace i by L 
in the conductivity formula, and hence get k ~ L. On the other hand for i >> L 
it is normally expected that a finite conductivity will be obtained. However, from 
the results presented in the previous sections, it is clear that this picture cannot 
be correct. At sufficiently large length scales, for low dimensional systems such as 
nanowires and nanotubes, we expect the conductivity to diverge as a power law 
k ~ L a . 

In the ballistic limit (defined as one where anharmonicity can be neglected) 
one can use the Landauer or NEGF formula and here there are examples where 
good agreement between theory and experiments can be seen. For nanotubes and 
nanowires with low impurity level it turns out that phonon mean free paths can 
be quite long and so transport is ballistic over fairly long length scales. 

Experiments on nanowires: One of the first measu reme nts of phonon thermal 



conductance of a nanosystem was that by Tighe et al. 1871 ] . In a beautiful exper- 
iment they measured the conductance of insulating GaAs wires of length ~ 5.5fim 
and cross-section ~ 200nmx 300nm. At low temperatures (1.5 — 5K) they obtained 
conductances of order ~ 10~ 9 W/K . From their data and using kinetic theory argu- 
ments they estimated the phonon mean free path to be of order ~ \[xm. Now one 
can ask the question: what is the thermal conductance of a perfectly transmitting 
ID wire? This can be easily obtained from Eq. (|65p by setting T(uj) = 1 and one 
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gets: 

G -Ar-2^y - ip ' (123) 

where we assume a phonon dispersion between — a; m . At temperatures T « 
hw/ks this gives: 

Gth = 50 = 7T 2 k 2 B T/(3h) (124) 

where g = (9.456 x 10~ 13 W/K*)T has been proposed as the quantum of ther- 
mal conductance 1881 ] and is the maximum value of energy transported per 



phonon mode. In a wire that is not strictly ID, as is the case of a wire of diameter 
~ 200nm, other modes would cont ribu te also to the current. In another nice and 
difficult experiment, Schwab et al. [l88( ] were able to measure the quantum of ther- 
mal conductance. They used a silicon nitride wire which had four lowest massless 
modes and other massive modes corresponding to the finite width of the wire. By 
going to sufficiently low temperatures (T < IK) they were able to suppress trans- 
port by the massive modes. Also one had to ensure very good contacts between 
the wires and reservoirs. The authors were able to verify that the resulting conduc- 
tance corresponded to the value go. The agreement is in fact quite impressive. At 
high temperatures all modes would contribute and so it will be difficult to verify 
the classical ID result with this system. Note that theoretically T(u) = 1 can 
probably be achieved only with Rubin baths. For other baths (for example Ohmic) 
this cannot be obtained even for ordered chains and as a result, the temperature 
dependence of conductivity can be quite different [for example see sec. (|3.3.ip ] from 
the linear dependence in E q. (|124[) . 



Another experiment 1891 ] reported measurements of thermal conductivity of sin- 
gle crystalline Si nanowires with wires several microns long and with varying di- 
ameters between 22nm to 115nm. They found that the thermal conductivity in- 
creased rapidly with diameter and was almost two orders of magnitude smaller 
(~ 20W/mK) than the bulk value. These results have not been clearly understood. 
So far there has been no experiments measuring the dependence of conductivity 
on length in nanowires. 

Experiments on nanotubes: Apart from nanowires, there have also been a 
number of measurements of heat current in nanotubes. One of the first measure- 
ments of conductance in individual samples was by Kim et al. fi~90l ]. on a 2Afim 
long and 14nm diameter multiwalled carbon nanotube (MWCNT). They found a 
very high thermal conductivity of ~ 3000W/mK (at room temperature) and noted 
that this corresponded to a phonon mean free of ~ 500nm. Somewh at su rprisingly 
this was close to a theoretically predicted value by Berber et al. [l9ll ] who had 
performed classical molecular dynamics simulations (using the Green-Kubo for- 
mula) for a (10, 10) carbon nanotube. Using realistic potentials they reported a 
high thermal conductivity of ~ 6000W/mK , at room temperature. From our ex- 
pectations of diverging conductivity we expect that these reported values, both in 
experiments and simulations, will i ncrease with increasing length of the wi re. I n 



fact simulations by Maruyama M, Zhang and Li & and by Yao et al. M 



again with realistic potentials, do find such a increase. A more recent simulation 
[196] however finds a converging conductivity. 
As far as theoretical work on heat conduction in ca rbon nanotubes is concerned, 



we mention the insightful paper by Mingo and Broido [1941 ] . who point out that it is 



necessary to perform quantum mechanical calculations in order to understand the 
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experimental results and that classical calculations can be misleading. They mostly 
consider the ballistic conductance using the Landauer formula [ Eq. (I65p ] and also 
the Boltzmann-Peierls equation at longer lengths. However they do not comment on 
the system size dependence of conductivity. It is likely that as far as the question of 
system size dependence of conductivity is concerned, the answer should probably be 
independent of whether one is doing a classical or a quantum calculation. Hence it 
will be useful to settle this issue through classical molecular dynamics simulations. 
From our experience with the difficulty in reaching asymptotic system sizes for 
ID and 2D systems, it is clear that one has to be careful before coming to quick 
conclusions. 



The experimental results of Fujii et al. [1971 ] on individual MWCNT give some 
hints of anomalous behaviour. They also obtain large values of thermal conductivity 
but find that it decreases with the diameter of the tubes. At room temperatures 
the thermal conductivity of a 3.7/um long nanotube with diameter ~ Wnm was 
about 2500W/mK while that of a 3.6/xm long, diameter ~ 30ram nanotube was 
about 500W/mK . Note that the dependence of k on diameter is opposite to that 
found for nanowires mentioned earlier. 

Measurements on indiv idua l single-wall carbon nanotubes (SWCNT) also have 



now been done. Yu et al. [1981 ] observed that the thermal conductance of a 2.76/xm 
long suspended tube was very close to the calculated ballistic thermal conductance 
(calculated using the Landauer formula) of a lnm diameter tube. In the tem- 
perature range of 100 — 300K they found increasing conductance and no signs of 
significant phonon-phonon sc atte ring. Another measurement on a single- walled car- 
bon nanotube by Pop et al. [199] measured the conductance to temperatures upto 
800K and they found a 1/T decay at large temperatures. They also report mea- 
surements on various lengths, ranging from 0.5/xm to 10/xm and diameter 1.5nm, 
and curiously, they found increasing conductance with length. This the authors ex- 
plain can be understood to be a result of the large phonon mean free path ~ 0.5/xm 
and phonon boundary scattering. 

An experimental proof of the divergence of thermal conductivity with system 
size is probably the dream of many theorists. There seems to be rapid progress in 
the direction of making this possible. The first indication of length dependence was 
reported by Wang et al. |200l ]. for the case of a SWCNT placed on a silicon sub- 
strate. They measured samples of lengths between 0.5 — 7fj,m at room temperature 
and found a slow increase of the conductance. The most recent experiments by 



Chang et al. [201 ] makes a detailed investigation of the length dependence of con- 



ductivity in multiwalled nanotubes, of carbon and boron-nitride, and claim to have 
found convincing evidence for violation of Fourier's law. These room temperature 
measurements were on suspended tubes of effective length between ~ 3.7 — 7/jm 
and their estimate of phonon mean free path is ~ 20 — 50nm. From their mea- 
surements (over the rather limited length scale) the authors conclude that a ~ 0.6 
for the carbon nanotube a ~ 0.5 for the boron nitride sample (which is isotopi- 
cally disordered). It is interesting to note that the two samples in this experiment, 
approximately correspond to the ordered and disordered FPU models. 
Experi men ts on suspended single layer graphene sheets have also been made 



recently 202] and so interesting experimental results from two dimensional systems 
can also be expected in the near future. It is of course too early to make definite 
conclusions from these experiments. 

Finally we briefly discuss one other area, that of thermal rectifiers, where an 
experiment was motivated by t heore tical work on simple models of heat conduction. 



In a paper by Terraneo et al. 2031 ]. an inhomogeneous nonlinear lattice model of 



heat conduction was proposed. This model had the interesting property that by 
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changing a single parameter on a part of the chain one could cause a transition from 
insulating to conducting behaviour. A related observation was that the absolute 
value of the heat current depended on the sign of AT = Tl — Tr. Thus one basically 
had a model for a thermal rectifier. Note that it can be proved rigorously, for both 
harmonic systems (with inhomogeneity but without self-consistent reservoirs) as 
well as homogeneous anharmonic systems, that J (AT) = — J (—AT) and so these 
systems cannot work as rectifiers. For harmonic systems this follows immediately 
from the general expression for current given in sec. (j3.2h . Clearly one needs both 
inhomogeneity as well as anharmonicity to get rectification and Terraneo et al. gave 
a simple explicit demonstration of how this could be achieved. Physically their 
results can be understood easily by thinking of the anharmonicity as giving rise 
to effective phonon bands which can be moved up and down by increasing or 
decreasing local temperatures. Phonon flow from the reservoirs into the system 
can thus be controll ed. 

Since the work in |203], a n umber of papers have observed t his effect in a number 
of models HqI, H S & Sol, [H, Eft [m|, HH, EH H • Based on the model 



of thermal recifier, Yang and Li have proposed a design of a thermal logic gate [2151 ] . 
An ex perimenta l observation of thermal rectification was made recently by Chang 
et al. |2ld . 217 ]. They made measurements of the heat current in a boron-nitride 



nanotube which was mass-loaded externally in an inhomogeneous way, and were 
able to obtain a small rectification. 

We conclude this section with the note that the situation looks hopeful for vig- 
orous interactions between theory and experiments. 



9. Concluding remarks 

The fact that Fourier's law is not valid in ID and 2D systems is a surprising result 
and probably the most important knowledge gained from the large number of stud- 
ies on heat conduction in low dimensional systems. Even in the limit of the system 
length being much larger than typical scattering lengths in a system, one finds that 
it is not possible to define a thermal conductivity as an intrinsic size-independent 
property of the system. This discovery is not only of academic interest but also 
important from the point of view of understanding real experiments. For example 
this tells us that it does not make sense to talk about the thermal conductivity of 
a carbon nanotube since this will keep changing with the length of the nanotube. 
To summarize, the main conclusions of this review are: 

(i) Fourier's law is not valid in momentum-conserving systems in one and two 
dimensions. 

For disordered harmonic systems, k, ~ N a , where a depends on boundary con- 
ditions and spectral properties of heat baths. 

For nonlinearly interacting systems without disorder, simulation results on a 
number of models indicate that in ID, a = 1/3, and that there is only one univer- 
sality class. There is disagreement between predictions from different theoretical 
approaches. In 2D, the theoretical prediction of n ~ log(A^) has not been verified 
in the latest simulations. 

(ii) Fourier's law, as far as the scaling J ~ 1/N is concerned, is valid in 
momentum-non-conserving non-integrable systems in all dimensions. Both theory 
and most simulations agree on this. 

(iii) In ID oscillator systems with both disorder and anharmonicity, the asymp- 
totic system size dependence of current is determined by anharmonicity alone, and 
localization becomes irrelevant. 

(iii) Chaos is neither a necessary nor a sufficient condition for validity of Fourier's 
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law. This result follows from the observation that Fourier's law is valid in billiard- 
like systems with polygonal scatterers which have zero-Lyapunov coefficient and 
hence are non-chaotic. On the other hand Fourier's law is not valid for the FPU 
system in any parameter regime even though it has positive Lyapunov exponents. 

(iv) For momentum-conserving ID systems, it is not possible to write Fourier's 
law in the form J = —KjyVT with kn defined as a size dependent conductivity. 
This follows from the anomalous steady state temperature profiles that seem to be 
invariably obtained in such systems. 

(v) For harmonic lattice systems, the Langevin equation Green's function 
(LEGF) formalism provides a very useful theoretical framework for understand- 
ing heat transport, in both classical and quantum systems. 

Some interesting problems that need to be addressed in the future are the fol- 
lowing: 

(a) Exact determination of the exponent a in any one dimensional momentum 
conserving model with purely Hamiltonian dynamics and without use of the usual 
Green-Kubo formula. 

(b) Simulations on nonlinearly interacting systems in ID, 2D and 3D for larger 
system sizes and different models, in order to establish the exponents convincingly. 
It will be nice to have more results on systems such as hard discs and spheres. 

We point out that the understanding of heat conduction even in three dimen- 
sional macroscopic syst ems is incomplete. A nice example of this can be seen from 
the discussion given in 2181 ] , in the context of understanding experimental results 



on heat conductivity of a highly purified single crystal diamond. A related point: 
in 3D it is a belief (see for example [93]) that at low temperatures, where Ump- 
klapp processes become exponentially rare, normal processes along with impurity 
scattering lead to a finite conductivity for the system. Can this be given some more 
rigorous justification, or, verified in simulations ? 

(c) Finding a for two and three dimensional disordered harmonic systems an- 
alytically. Further simulations are also necessary here. What are the connections 
with localization theory ? 

(d) For disordered anharmonic systems, for disorder strength and anharmonicity 
strengths denoted by A and A respectively, what is the phase diagram in the A — A 
plane ? 

(e) For open systems there is a rigorous derivation of a linear response result 
which is valid for finite systems. Is it possible to prove the equivalence of this 
with the usual Green-Kubo formula for closed systems, in some example ? This is 
probably true for systems with normal transport and probably not true for systems 
with anomalous transport. 

(f) Proof of non-existence (or existence) of phase transitions in one dimensional 
models of heat conduction with short range interactions. 

(g) One needs studies for quantum interacting systems since most experimental 
work seems to be in this domain. We have seen that the Green's function approach 
has been successful in understanding harmonic systems (i.e. ballistic transport). 
An extension of this approach to the anharmonic case would be very useful. Apart 
from the Green's function formalism, the approach used by Chen et al. [1^] may 
be a useful method for this problem. 

(h) How valid is the hydrodynamic description for systems with anomalous trans- 
port? If they are valid, what are the correct hydrodynamic equations? For example 
we have seen that one cannot use the equation J = —kn(T)VT to describe the 
steady state. 

(i) Non steady state properties: most of the studies on heat conduction have been 
on measurement of current and temperature in the nonequilibrium steady state. 



November 19, 2008 19:56 Advances in Physics reva 



72 REFERENCES 

In general of course one is interested also in time dependent properties. In fact the 
diffusion equation following from Fourier's law is itself a time-dependent equation. 
Also many experiments make measurements in non-steady state conditions, such 
as by studying heat pulses and frequency dependent studies. Thus it is necessary to 
have more theoretical studies on heat flow in non steady state situations. Interesting 
questions can be asked here, e.g., can one talk of a frequency dependent thermal 



conductivity [219 ] 



? 

(j) Scalar versus vector models: for lattice models one question is whether the 
dimensionality of the displacement vectors matters as far as exponents are con- 
cerned. While it is usually assumed that dimensionality does not matter, it will be 
nice to have a proof of this. 

(k) Temperature dependence of conductivity or conductance: Apart from the 
system size dependence it will be useful to get more results on the temperature 
dependence of the linear response heat current, since this is one of the things that 
experimentalists are interested in. 
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